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

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