Загрузите модель роста опухоли.

sbioloadproject tumor_growth_vpop_sa.sbproj

Заставьте вариант предполагаемыми параметрами и дозой применяться к модели.

v = getvariant(m1);
d = getdose(m1,'interval_dose');

Получите активный configset и установите вес опухоли как ответ.

cs = getconfigset(m1);
cs.RuntimeOptions.StatesToLog = 'tumor_weight';

Симулируйте модель и постройте профиль роста опухоли.

sbioplot(sbiosimulate(m1,cs,v,d));

Выполните глобальный анализ чувствительности (GSA) на модели, чтобы найти параметры модели, к которым рост опухоли чувствителен.

Во-первых, получите параметры модели интереса, которые вовлечены в фармакодинамику роста опухоли. Задайте ответ модели как вес опухоли.

modelParamNames = {'L0','L1','w0','k1','k2'};
outputName = 'tumor_weight';

Затем выполните GSA путем вычисления первого - и общий порядок индексы Sobol с помощью sbiosobol. Установите 'ShowWaitBar' к true показать прогресс симуляции. По умолчанию функция использует 1 000 выборок параметра, чтобы вычислить индексы Sobol [1].

sobolResults = sbiosobol(m1,modelParamNames,outputName,'Variants',v,'Doses',d,'ShowWaitBar',true)

sobolResults = 
  Sobol with properties:

                Time: [444x1 double]
        SobolIndices: [5x1 struct]
            Variance: [444x1 table]
         Observables: {'[Tumor Growth Model].tumor_weight'}
    ParameterSamples: [1000x5 table]
      SimulationInfo: [1x1 struct]

Можно изменить количество выборок путем определения 'NumberSamples' аргумент пары "имя-значение". Функция требует в общей сложности (number of input parameters + 2) * NumberSamples симуляции модели.

Покажите средний ответ модели, результаты симуляции и теневую область, покрывающую 90% результатов симуляции.

plotData(sobolResults);

Можно настроить область квантиля к различному проценту путем определения 'Alphas' для более низких и верхних квантилей всех ответов модели. Например, альфа-значение 0,1 графиков теневая область между 100 * alpha и 100 * (1 - alpha) квантили всех симулированных ответов модели.

plotData(sobolResults,'Alphas',0.1);

Постройте курс времени первого - и общий порядок индексы Sobol.

h = plot(sobolResults);
% Resize the figure.
h.Position(:) = [100 100 1280 800];

Индекс Sobol первого порядка входного параметра дает часть полного отклонения ответа, которое может быть приписано изменениям одного только входного параметра. Индекс общего порядка дает часть полного отклонения ответа, которое может быть приписано любым объединенным изменениям параметра, которые включают изменения входного параметра.

Из графиков индексов Sobol, параметры L1 и w0 кажется, самые чувствительные параметры к весу опухоли, прежде чем доза была применена в t = 7. Но после того, как доза применяется, k1 и k2 станьте более чувствительными параметрами и способствуйте больше всего дозируемому после того, как этапу веса опухоли. Общий график отклонения также показывает большее отклонение для этапа после-того,-как-дозы в t> 35, чем для этапа прежде-чем-дозы роста опухоли, указывая на тот k1 и k2 могут быть более важные параметры, чтобы заняться расследованиями далее. Часть необъясненного отклонения показывает некоторое отклонение приблизительно в t = 33, но общий график отклонения показывает мало отклонения в t = 33, означая, что необъясненное отклонение могло быть незначительным. Часть необъясненного отклонения вычисляется как 1 - (сумма всех индексов Sobol первого порядка), и общее отклонение вычисляется с помощью var(response), где response ответ модели в каждом моменте времени.

Можно также отобразить величины чувствительности в столбиковой диаграмме.

bar(sobolResults)

Можно задать больше выборок, чтобы увеличить точность индексов Sobol, но симуляция может занять больше времени, чтобы закончиться. Используйте addsamples добавить больше выборок. Например, если вы задаете 1 500 выборок, функция выполняет 1500 * (2 + number of input parameters) симуляции.

gsaMoreSamples = addsamples(gsaResults,1500)

Свойство SimulationInfo объекта результата содержит различную информацию для вычисления индексов Sobol. Например, данные о симуляции модели (SimData) для каждой симуляции с помощью набора выборок параметра хранятся в SimData поле свойства. Это поле является массивом SimData объекты.

sobolResults.SimulationInfo.SimData

 
   SimBiology SimData Array : 1000-by-7
 
   Index:    Name:         ModelName:         DataCount: 
   1           -           Tumor Growth Model 1          
   2           -           Tumor Growth Model 1          
   3           -           Tumor Growth Model 1          
   ...                                                   
   7000        -           Tumor Growth Model 1          
 

Можно узнать, перестала ли какая-либо симуляция модели работать во время расчета путем проверки ValidSample поле SimulationInfo. В этом примере поле не показывает не пройдено запусков симуляции.

all(sobolResults.SimulationInfo.ValidSample)

ans = 1x7 logical array

   1   1   1   1   1   1   1

SimulationInfo.ValidSample таблица логических значений. Это имеет тот же размер как SimulationInfo.SimData. Если ValidSample указывает, что любые симуляции перестали работать, можно получить больше информации о тех запусках симуляции и выборках, используемых для тех запусков путем извлечения информации из соответствующего столбца SimulationInfo.SimData. Предположим, что четвертый столбец содержит один или несколько не пройдено запусков симуляции. Получите данные моделирования и демонстрационные значения, используемые для той симуляции с помощью getSimulationResults.

[samplesUsed,sd,validruns] = getSimulationResults(sobolResults,4);

Можно добавить пользовательские выражения как observables и вычислить индексы Sobol для добавленного observables. Например, можно вычислить индексы Sobol для максимального веса опухоли путем определения пользовательского выражения можно следующим образом.

% Suppress an information warning that is issued during simulation.
warnSettings = warning('off', 'SimBiology:sbservices:SB_DIMANALYSISNOTDONE_MATLABFCN_UCON');
% Add the observable expression.
sobolObs = addobservable(sobolResults,'Maximum tumor_weight','max(tumor_weight)','Units','gram');

Постройте вычисленные результаты симуляции, показывающие 90%-ю область квантиля.

h2 = plotData(sobolObs);
h2.Position(:) = [100 100 1280 800];

Можно также удалить заметное путем определения его имени.

gsaNoObs = removeobservable(sobolObs,'Maximum tumor_weight');

Восстановите настройки предупреждения.

warning(warnSettings);

Загрузите модель Установленного целью расположения препарата (TMDD).

sbioloadproject tmdd_with_TO.sbproj

Получите активный configset и установите целевое заполнение (TO) как ответ.

cs = getconfigset(m1);
cs.RuntimeOptions.StatesToLog = 'TO';

Симулируйте модель и постройте TO профиль.

sbioplot(sbiosimulate(m1,cs));

Задайте воздействие (область под кривой профиля TO) порог для целевого заполнения.

classifier = 'trapz(time,TO) <= 0.1';

Выполните MPGSA, чтобы найти чувствительные параметры относительно TO. Варьируйтесь значения параметров между предопределенными границами, чтобы сгенерировать 10 000 выборок параметра.

% Suppress an information warning that is issued during simulation.
warnSettings = warning('off', 'SimBiology:sbservices:SB_DIMANALYSISNOTDONE_MATLABFCN_UCON');
rng(0,'twister'); % For reproducibility
params = {'kel','ksyn','kdeg','km'};
bounds = [0.1, 1; 
          0.1, 1;
          0.1, 1;
          0.1, 1];
mpgsaResults = sbiompgsa(m1,params,classifier,'Bounds',bounds,'NumberSamples',10000)

mpgsaResults = 
  MPGSA with properties:

                    Classifiers: {'trapz(time,TO) <= 0.1'}
    KolmogorovSmirnovStatistics: [4x1 table]
                       ECDFData: {4x4 cell}
              SignificanceLevel: 0.0500
                        PValues: [4x1 table]
              SupportHypothesis: [10000x1 table]
                    Observables: {'TO'}
               ParameterSamples: [10000x4 table]
                 SimulationInfo: [1x1 struct]

Постройте квантили симулированного ответа модели.

plotData(mpgsaResults);

Постройте эмпирические кумулятивные функции распределения (eCDFs) принятых и отклоненных выборок. За исключением km, ни один из параметров не показывает значительную разницу в eCDFs для принятых и отклоненных выборок. km постройте показывает крупному Кольмогорову-Смирнову (K-S) расстояние между eCDFs принятых и отклоненных выборок. Расстояние K-S является максимальным абсолютным расстоянием между двумя кривыми eCDFs.

h = plot(mpgsaResults);
% Resize the figure.
pos = h.Position(:);
h.Position(:) = [pos(1) pos(2) pos(3)*2 pos(4)*2];

Чтобы вычислить расстояние K-S между двумя eCDFs, SimBiology использует двухсторонний тест на основе нулевой гипотезы, что два распределения принятых и отклоненных выборок равны. Смотрите kstest2 (Statistics and Machine Learning Toolbox) для деталей. Если расстояние K-S является большим, то эти два распределения отличаются, означая, что классификация выборок чувствительна к изменениям входного параметра. С другой стороны, если расстояние K-S мало, то изменения входного параметра не влияют на классификацию выборок. Результаты предполагают, что классификация нечувствительна к входному параметру. Чтобы оценить значение статистической величины K-S, отклоняющей нулевую гипотезу, можно исследовать p-значения.

bar(mpgsaResults)

Столбиковая диаграмма показывает две панели для каждого параметра: один для расстояния K-S (статистическая величина K-S) и другой для соответствующего p-значения. Вы отклоняете нулевую гипотезу, если p-значение меньше уровня значения. Крест (x) показан для любого p-значения, которое является почти 0. Вы видите, что точное p-значение соответствует каждому параметру.

[mpgsaResults.ParameterSamples.Properties.VariableNames',mpgsaResults.PValues]

ans=4×2 table
      Var1      trapz(time,TO) <= 0.1
    ________    _____________________

    {'kel' }          0.0021877      
    {'ksyn'}                  1      
    {'kdeg'}            0.99983      
    {'km'  }                  0      

P-значения km и kel меньше уровня (0.05) значения, поддерживая альтернативную гипотезу, что принятые и отклоненные выборки прибывают из различных распределений. Другими словами, классификация выборок чувствительна к km и kel но не к другим параметрам (kdeg и ksyn).

Можно также построить гистограммы принятых и отклоненных выборок. historgrams позволяют вам видеть тренды в принятых и отклоненных выборках. В этом примере, гистограмме km показывает, что существуют более принятые выборки для большего km значения, в то время как kel гистограмма показывает, что существует меньше отклоненных выборок как kel увеличения.

h2 = histogram(mpgsaResults);
% Resize the figure.
pos = h2.Position(:);
h2.Position(:) = [pos(1) pos(2) pos(3)*2 pos(4)*2];

Восстановите настройки предупреждения.

warning(warnSettings);