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

sbioloadproject tumor_growth_vpop_sa.sbproj

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

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

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

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, вычислив индексы Соболя первого и общего порядков с помощью sbiosobol. Задайте 'ShowWaitBar' на true чтобы показать прогресс симуляции. По умолчанию функция использует 1000 выборок параметров, чтобы вычислить индексы Соболь [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);

Постройте график хода индексов Соболь первого и общего порядков.

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

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

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

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

bar(sobolResults)

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

gsaMoreSamples = addsamples(gsaResults,1500)

Свойство SimulationInfo объекта результата содержит различную информацию для вычисления индексов Соболь. Например, данные симуляции модели (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.SimDatА. Предположим, что четвертый столбец содержит один или несколько неудачных запуски симуляции. Получите данные моделирования и значения выборки, используемые для этой симуляции, используя getSimulationResults.

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

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

% 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

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

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

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

sbioplot(sbiosimulate(m1,cs));

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

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

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

% 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ни один из параметров не показывает значительного различия в eCDF для принятых и отклонённых выборок. The km график показывает большое расстояние Колмогоров-Смирнов (К-С) между eCDF принятых и отклонённых выборок. Расстояние 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 между двумя eCDF, 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).

Можно также построить гистограммы принятых и отклонённых выборок. Исторические граммы позволяют вам увидеть тренды в принятых и отклонённых выборках. В этом примере гистограмма 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);