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

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];

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

Из индексов Соболь сюжеты, параметры 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)

Свойство ReservateInfo объекта результата содержит различную информацию для вычисления индексов Соболь. Например, данные моделирования модели (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 для поиска чувствительных параметров относительно транспортного заказа. Изменение значений параметров между предварительно определенными границами для создания 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);

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

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 (Статистика и инструментарий машинного обучения) для получения подробной информации. Если расстояние 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);