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

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