Ссылка

Количественная характеристика дрожжей гетеротримерный цикл белка G. Тау-Му И, Ироаки Китано и Мелвин Ай. Саймон. PNAS (2003) издание 100, 10764-10769.

Цели

Фон

В дрожжах Saccharomyces cerevisiae, G белок, сигнализирующий в сцепляющемся ответе, является хорошо охарактеризованной трассой преобразования сигнала. Феромон, спрятавший альфа-ячейками, активируется, G-белок связал факторный альфой приемник (Ste2p) в ячейки, который приводит ко множеству клеточных реакций включая фиксацию клеточного цикла и синтез новых белков. G белки и белок G связался, приемники (GPCRs) являются особым вниманием усилий по изобретению лекарства в фармацевтической промышленности. Много проданных наркотиков предназначаются для GPCRs - некоторые примеры включают, те для сокращения кислоты желудочного сока (ранитидин, целевой гистаминовый приемник H2), мигрень (sumatriptan, предназначается для подтипа приемника серотонина), шизофрения (olanzapine, целевой серотонин и приемники допамина), и аллергии (дезлоратадин, целевые гистаминовые приемники). Далее, некоторые оценки предполагают, что GPCRs являются целенаправленным особым вниманием 40% усилий по изобретению лекарства. Один подход является к модели GPCR сигнальными трассами, чтобы анализировать и предсказать и нисходящие эффекты и эффекты в связанных трассах. Этот пример исследует построение моделей, симуляцию и анализ цикла белка G в трассе ответа феромона дрожжей.

Этот рисунок является графическим представлением концептуальной основы, используемой, чтобы смоделировать дрожжи G цикл белка.

Следующие сокращения используются:

Реакции трассы

Как показано на рисунке, цикл может быть сжат в набор биохимических реакций:

1) Взаимодействие лиганда приемника (обратимая реакция)

  L + R < - > RL

2) Гетеротримерное формирование G-белка

  Gd + Gbg -> G

3) Активация G-белка - Примечание, ниже которого RL появляется с обеих сторон уравнения, потому что RL является модификатором или катализатором для реакции. RL ни не производится, ни используется этой реакцией.

  RL + G -> Ga + Gbg + RL

4) Синтез приемника и ухудшение (обработанный как обратимая реакция представлять ухудшение и синтез)

  R < - > null

5) Ухудшение лиганда приемника

  RL -> null

6) Деактивация G-белка, катализируемая Sst2p в диком типе, напрягает TMY101, и некатализируемый в мутанте напрягают TMY111 с разрушением в гене SST2.

  Ga -> Gd

Все значения были преобразованы в молекулу для сумм разновидностей и молекулу/секунда, или 1/второй для параметров уровня.

Создавание модели SimBiology® для трассы Дикого Типа

Создайте объект модели SimBiology с именем 'Гетеротримерный вес Белка G'.

 modelObj = sbiomodel('Heterotrimeric G Protein wt');

Добавьте взаимодействие Лиганда Приемника (обратимая реакция).

 reactionObj1 = addreaction(modelObj, 'L + R <-> RL', ...
     'Name', 'Receptor-ligand interaction');

Используйте 'MassAction' кинетический закон для реакции. Эта модель создана с помощью кинетики действующих масс для всех реакций.

 kineticlawObj1 = addkineticlaw(reactionObj1, 'MassAction');

Добавьте прямые и противоположные параметры уровня.

 addparameter(modelObj, 'kRL', 3.32e-18);
 addparameter(modelObj, 'kRLm', 0.01);

Присвойте ParameterVariableNames в кинетическом объекте закона. Это сопоставляет ParameterVariables с ParameterVariableNames в кинетическом объекте закона так, чтобы скорость реакции могла быть определена.

kineticlawObj1.ParameterVariableNames = {'kRL', 'kRLm'};

SimBiology автоматически создает объекты разновидностей для каждой из участвующих разновидностей в реакциях. Установите начальные суммы этих разновидностей.

% Set initial amount for 'L'
 modelObj.Reactions(1).Reactants(1).InitialAmount = 6.022E17;
% Set initial amount for 'R'
 modelObj.Reactions(1).Reactants(2).InitialAmount = 10000.0;
% Leave initial amount for 'RL' at default value (0.0)

ReactionRate для первой реакции был теперь сконфигурирован.

 reactionObj1.ReactionRate

ans =

    'kRL*L*R - kRLm*RL'

Завершение модели Дикого Типа

Чтобы создать модель деформации дикого типа (TMY101), добавьте остальную часть реакций и параметров, создайте кинетические объекты закона для каждой из реакций и присвойте переменные параметра для кинетических законов.

Добавьте и сконфигурируйте реакцию для Гетеротримерного формирования G-белка.

 reactionObj2 = addreaction(modelObj, 'Gd + Gbg -> G', ...
     'Name', 'G protein complex formation');
 kineticlawObj2 = addkineticlaw(reactionObj2, 'MassAction');
 addparameter(modelObj, 'kG1', 1.0);
 kineticlawObj2.ParameterVariableNames = 'kG1';
% Set initial amount for 'Gd'
 modelObj.Reactions(2).Reactants(1).InitialAmount = 3000;
% Set initial amount for 'Gbg'
 modelObj.Reactions(2).Reactants(2).InitialAmount = 3000;
% Set initial amount for 'G'
 modelObj.Reactions(2).Products(1).InitialAmount = 7000;

Добавьте и сконфигурируйте реакцию для активации G-белка.

 reactionObj3 = addreaction(modelObj, 'G + RL -> Ga + Gbg + RL', ...
     'Name', 'G protein activation');
 kineticlawObj3 = addkineticlaw(reactionObj3, 'MassAction');
 addparameter(modelObj, 'kGa', 1.0E-5);
 kineticlawObj3.ParameterVariableNames = 'kGa';
% Set initial amount for 'Ga'
 modelObj.Reactions(3).Products(1).InitialAmount =  0.0;

Добавьте и сконфигурируйте реакцию для синтеза приемника и ухудшения.

 reactionObj4 = addreaction(modelObj, 'R <-> null', ...
     'Name', 'R synthesis/degradation');
 kineticlawObj4 = addkineticlaw(reactionObj4, 'MassAction');
 addparameter(modelObj, 'kRdo', 4.0E-4);
 addparameter(modelObj, 'kRs', 4.0);
 kineticlawObj4.ParameterVariableNames = {'kRdo','kRs'};

Добавьте и сконфигурируйте реакцию для ухудшения лиганда приемника.

 reactionObj5 = addreaction(modelObj, 'RL -> null', 'Name', 'RL degradation');
 kineticlawObj5 = addkineticlaw(reactionObj5, 'MassAction');
 addparameter(modelObj, 'kRD1', 0.0040);
 kineticlawObj5.ParameterVariableNames = 'kRD1';

Добавьте и сконфигурируйте реакцию для деактивации G-белка.

 reactionObj6 = addreaction(modelObj, 'Ga -> Gd', 'Name', 'Gprotein inactivation');
 kineticlawObj6 = addkineticlaw(reactionObj6, 'MassAction');
 addparameter(modelObj, 'kGd', 0.11);
 kineticlawObj6.ParameterVariableNames = 'kGd';

Проверяйте ReactionRate всех реакций.

 get(modelObj.Reactions, {'Reaction', 'ReactionRate'})

ans =

  6x2 cell array

    {'L + R <-> RL'           }    {'kRL*L*R - kRLm*RL'}
    {'Gd + Gbg -> G'          }    {'kG1*Gd*Gbg'       }
    {'G + RL -> Ga + Gbg + RL'}    {'kGa*G*RL'         }
    {'R <-> null'             }    {'kRdo*R - kRs'     }
    {'RL -> null'             }    {'kRD1*RL'          }
    {'Ga -> Gd'               }    {'kGd*Ga'           }

Симуляция модели Дикого Типа и графический вывод результатов

Чтобы отметить быстрое повышение и последующее снижение разновидностей Ga, симулируйте модель в течение 600 с и сохраните результаты.

Измените StopTime объекта конфигурации модели по умолчанию с 10-х (simulationTime) к 600 с. Кроме того, не делайте данных логов для лиганда 'L' (modelObj.Species (1)), потому что это берет значения, которые являются порядками величины выше, чем другие разновидности. Это делает визуализацию разновидностей в графике более удобной. Чтобы выполнить это, задайте StatesToLog, чтобы включать все разновидности кроме 'L'.

 configsetObj = getconfigset(modelObj);
 configsetObj.StopTime = 600;
 configsetObj.SolverOptions.AbsoluteTolerance = 1.e-9;
 configsetObj.RuntimeOptions.StatesToLog = ...
     sbioselect(modelObj, 'Type', 'species', 'Where', 'Name', '~=', 'L');

Симулируйте модель и возвратите результаты к этим трем переменным 'время', 'данные' и 'имена'.

 [time, data, names] = sbiosimulate(modelObj);

Отобразите данные на графике.

 plot(time, data);
 legend(names, 'Location', 'NorthEastOutside');
 xlabel('Time (seconds)');
 ylabel('Species Amounts');
 grid on;

Создание варианта модели для деформации мутанта

Модель цикла G-белка для деформации мутанта отличается по уровню, на котором происходит деактивация активного G-белка (Ga). Этим уровнем управляет значение параметра уровня kGd. Можно представлять деформацию мутанта с помощью объекта Variant. Вариант SimBiology хранит альтернативные значения для одного или нескольких свойств модели SimBiology, таких как InitialAmount разновидности или Значение параметра.

Добавьте вариант, названный 'мутантом' к модели.

 variantObj = addvariant(modelObj, 'mutant');

Добавьте содержимое в вариант, чтобы задать альтернативное значение 0,004 для параметра kGd.

 addcontent(variantObj, {'parameter', 'kGd', 'Value', 0.004});

Симуляция трассы мутанта и графический вывод результатов

Симулируйте модель с помощью объекта варианта мутанта. Это применяет значение 0,004 к параметру kGd в процессе моделирования. Возвратите результаты симуляции в объекте SimData. В дополнение к хранению данных моделирования SimBiology объекты SimData предоставляют методы для доступа к данным, графического вывода и анализа.

Установите свойство Active объекта варианта мутанта к истине и симулируйте.

 variantObj.Active = true;
 mutantData = sbiosimulate(modelObj);

Отобразите данные на графике с помощью пунктирных линий. См. также sbioplot для удобного графического вывода объектов SimData.

 plot(mutantData.Time, mutantData.Data, 'LineStyle', '--');
 legend(mutantData.DataNames, 'Location', 'NorthEastOutside');
 xlabel('Time (seconds)');
 ylabel('Species Amounts');
 grid on;

Сравните поведение активных разновидностей G-Protein (Ga) в трассах мутанта и диком типе.

 GaIndex = strcmp(names, 'Ga'); % index for wild-type results
 [tmut, xmut] = selectbyname(mutantData, 'Ga');
 plot(time, data(:,GaIndex), tmut, xmut, '--');
 xlabel('Time (seconds)');
 ylabel('Species Amounts');
 legend({'Ga (wt)','Ga (mutant)'}, 'Location', 'NorthEastOutside');
 grid on;

Выполнение скана параметра

Уровень деактивации G-белка намного ниже в деформации мутанта относительно дикого типа (kGd = 0.004 по сравнению с kGd = 0.11), который объясняет более высокие уровни активированного G-белка (Ga), в зависимости от времени наблюдаемый в вышеупомянутом сравнении. Для более подробного взгляда на то, как изменение kGd влияет на уровни Ga, выполните скан параметра нескольких симуляций, по которым значение kGd варьируется в области значений значений. Следующий пример иллюстрирует, что параметр сканирует более чем пять значений kGd; чтобы увеличить число итераций, измените значения в аргументах для linspace функция ниже.

Сгенерируйте пять равномерно распределенных kGd значений в пределах от 0,001 к 0,15.

 kGdValues = linspace(1e-3, 0.15, 5);

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

 scanData = [];

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

 sbioaccelerate(modelObj);

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

for kGd = kGdValues
    % Set the desired value of kGd in the variant.
    variantObj.Content{1}{4} = kGd;

    % Simulate the model, storing the results in a SimData object.
    sd = sbiosimulate(modelObj);
    scanData = [scanData; sd];
end

scanData теперь пять массивов элемента объектов SimData. Каждый объект содержит данные из одного запуска в скане параметра.

Извлеките timecourses для Ga от объектного массива SimData и графика на одной оси. Следующий код создает постепенный график; альтернативно, смотрите sbioplot и sbiosubplot.

 [tscan, xscan] = selectbyname(scanData, 'Ga');

 fh = figure;
 hold on;
 for c = 1:5
    plot(tscan{c}, xscan{c});
    str = sprintf(' k = %5.3f', kGdValues(c));
    text(tscan{c}(end), xscan{c}(end), str);
 end

 % Annotate the plot.
 axis(gca(fh), 'square');
 title('Varying the Value of kGd: Effect on Ga');
 xlabel('Time (seconds)');
 ylabel('Species Amounts');
 grid on;
 hold off;

Оценка параметра - фон

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

В этом разделе примера мы исследуем функциональность оценки параметра в контексте попытки подбирать модель белка G к экспериментальным данным.

Оценка параметра - сравнение результатов модели к экспериментальным данным

Для экспериментальных данных Рис. 5 ссылочной бумаги содержит timecourse для части активного белка G.

Сохраните экспериментальное время и утвердите данные.

  tExpt = [0 10 30 60 110 210 300 450 600]';
  GaFracExpt = [0 0.35 0.4 0.36 0.39 0.33 0.24 0.17 0.2]';
  data = groupedData(table(tExpt, GaFracExpt));
  data.Properties.IndependentVariableName = 'tExpt';

Вместо того, чтобы преобразовать это экспериментальные данные в абсолютные количества Ga, добавьте эту часть в модель с помощью непостоянного параметра и правила repeatedAssignment.

 GaFracObj = modelObj.addparameter('GaFrac', 'ConstantValue', 0);
 GaFracRule = modelObj.addrule('GaFrac = Ga / (Ga + G + Gd)', 'repeatedAssignment')

GaFracRule = 

   SimBiology Rule Array

   Index:    RuleType:             Rule:                      
   1         repeatedAssignment    GaFrac = Ga / (Ga + G + Gd)

Измените RuntimeOptions на конфигурации модели, чтобы регистрировать GaFrac.

 configsetObj.RuntimeOptions.StatesToLog = GaFracObj;

Деактивируйте вариант мутанта.

 variantObj.Active = false;

Симулируйте модель, храня результаты в объекте SimData.

 sdWild = sbiosimulate(modelObj);

Заставьте данные для 'GaFrac' использоваться позже в графике.

 [tWild, GaFracWild] = selectbyname(sdWild, 'GaFrac');

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

 sdWildResampled = resample(sdWild, tExpt, 'pchip');

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

 [~, GaFracWildResampled] = selectbyname(sdWildResampled, 'GaFrac');

Вычислите значение R-квадрата, измеряющее подгонку между симулированными и экспериментальными данными.

 sst = norm(GaFracExpt - mean(GaFracExpt))^2;
 sse = norm(GaFracExpt - GaFracWildResampled)^2;
 rSquare = 1-sse/sst;

Постройте результаты симуляции против экспериментальных данных для Ga.

 fh = figure;
 plot(tExpt, GaFracExpt, 'ro');
 legendText = {'Experiment'};
 title('Fit to Experimental Data for GaFrac');
 xlabel('Time (seconds)');
 ylabel('Species Amount');
 hold on;
 plot(tWild, GaFracWild);
 legendText{end+1} = sprintf('Original, R^2 = %4.2f', rSquare);
 legend(legendText{:});
 grid on;

Оценка параметра - оценка одного параметра модели

От скана параметра мы видели, что значение параметра kGd оказывает значительное влияние на timecourse разновидностей Ga. Давайте смотреть, можем ли мы улучшить подгонку результатов модели к экспериментальным данным путем варьирования значения kGd.

Выполните оценку параметра против экспериментальных данных, оптимизировав значение kGd. Постройте информацию об итерациях, в то время как оптимизация прогрессирует максимум до 15 итераций.

 paramToEst = estimatedInfo('kGd');
 kGdObj = sbioselect(modelObj, 'Name', 'kGd');
 opt = optimset('PlotFcns',@optimplotfval,'MaxIter',15);
 result1 = sbiofit(modelObj, data, 'GaFrac = GaFracExpt', paramToEst, ...
     [], 'fminsearch', opt);
 estValues1 = result1.ParameterEstimates

estValues1 =

  1x3 table

     Name      Estimate    StandardError
    _______    ________    _____________

    {'kGd'}    0.12142       0.0018548  

Сохраните ориентировочную стоимость kGd в новой модели Variant.

 optimVariantObj = addvariant(modelObj, 'Optimized kGd');
 addcontent(optimVariantObj, {'parameter', 'kGd', 'Value', estValues1.Estimate});

Активируйте новый вариант и деактивируйте вариант 'мутанта'.

 optimVariantObj.Active = true;
 mutantVariantObj = getvariant(modelObj, 'mutant');
 mutantVariantObj.Active = false;

Симулируйте модель с помощью ориентировочной стоимости kGd.

 sdEst1 = sbiosimulate(modelObj);

Отобразите данные на графике для GaFrac и сравните с предыдущими результатами.

 [t1, GaFracEst1] = selectbyname(sdEst1, 'GaFrac');
 sdEst1Resampled = resample(sdEst1, tExpt, 'pchip');
 [~, GaFracEst1Resampled] = selectbyname(sdEst1Resampled, 'GaFrac');
 sse1 = norm(GaFracExpt - GaFracEst1Resampled)^2;
 rSquare1 = 1-sse1/sst;
 figure(fh);
 plot(t1, GaFracEst1, 'm-');
 legendText{end+1} = sprintf('kGd Changed, R^2 = %4.2f', rSquare1);
 legend(legendText{:});

От значений R-квадрата мы видим, что подгонка к экспериментальным данным немного лучше с новой, ориентировочной стоимостью kGd. Если бы исходное значение для kGd было только грубой оценкой, мы могли бы интерпретировать эти результаты или как подтверждение первоначальной оценки или как улучшение по сравнению с ним.

Анализ чувствительности - фон

До сих пор мы интересовались динамическим поведением активного белка G, разновидностей Ga. Скан параметра показал, что эта разновидность значительно затронута значением константы скорости kGd управляющая деактивация G-белка. Используя оценку параметра, мы нашли, что путем оптимизации значения kGd, смогли лучше соответствовать экспериментальному timecourse для Ga.

Естественный вопрос спросить, что другие параметры влияния модели уровни Ga, и каковы величины тех эффектов? Анализ чувствительности позволяет вам отвечать на эти вопросы путем вычисления зависящих от времени производных одной или нескольких разновидностей ("выходные параметры") или относительно значений параметра модели или относительно начальных условий разновидностей ("входные факторы").

Анализ чувствительности - вычисление чувствительности

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

Деактивируйте объект варианта мутанта на модели так, чтобы чувствительность была вычислена с kGd в ее исходном значении.

 optimVariantObj.Active = false;

Настройте вычисление чувствительности в модели configset.

% Turn on SensitivityAnalysis in the solver options.
 configsetObj.SolverOptions.SensitivityAnalysis = true;

% Configure the sensitivity outputs and inputs for sensitivity analysis.
 sensitivityOpt = configsetObj.SensitivityAnalysisOptions;
 GaObj = sbioselect(modelObj, 'Type', 'species', 'Name', 'Ga');
 sensitivityOpt.Outputs = GaObj;
 params = sbioselect(modelObj, 'Type', 'parameter', 'Where', 'Name', '~=', 'GaFrac');
 sensitivityOpt.Inputs = params;
 sensitivityOpt.Normalization = 'Full';

Симулируйте модель.

 sdSens = sbiosimulate(modelObj);

Извлеките данные о чувствительности из объекта SimData и постройте вычисленную чувствительность.

 [t, R, sensOutputs, sensInputs] = getsensmatrix(sdSens);
 R = squeeze(R);

 figure;
 plot(t,R);
 title('Normalized sensitivity of Ga with respect to various parameters');
 xlabel('Time (seconds)');
 ylabel('Sensitivity');
 legend(sensInputs, 'Location', 'NorthEastOutside');
 grid on;

Оценка параметра - оценка нескольких параметров

Эти результаты показывают, что Ga не только чувствителен к параметру kGd, но также и к kGa, kRs, и kRD1. (Другая чувствительность неотличима от нуля на графике.) Варьирующийся эти четыре параметра могут сделать подгонку к экспериментальным данным лучше все еще.

Оцените, что эти четыре параметра совпадают с целевыми данными. Используйте ранее сконфигурированные опции оптимизации и текущие значения параметров в модели как начальная точка для оптимизации.

Выберите параметры kGa, kRs, kRD1, и kGd для оценки.

 paramsToEst = estimatedInfo({'kGa', 'kRs', 'kRD1', 'kGd'});

Оценка параметра проигнорирует опцию анализа чувствительности, если это будет включено в configset. Выключите SensitivityAnalysis в опциях решателя, чтобы избежать предупреждений.

 configsetObj.SolverOptions.SensitivityAnalysis = false;
 result2 = sbiofit(modelObj, data, 'GaFrac = GaFracExpt', paramsToEst, ...
     [], 'fminsearch', opt);
 estValues2 = result2.ParameterEstimates

estValues2 =

  4x3 table

      Name       Estimate     StandardError
    ________    __________    _____________

    {'kGa' }    9.0068e-06     3.0249e-06  
    {'kRs' }         4.549         11.786  
    {'kRD1'}     0.0031018      0.0027417  
    {'kGd' }       0.12381       0.053702  

Сохраните ориентировочные стоимости этих четырех параметров в новом варианте модели.

 optimVariantObj2 = addvariant(modelObj, 'Four parameter optimization');
 addcontent(optimVariantObj2, {'parameter','kGa', 'Value', estValues2.Estimate(1)});
 addcontent(optimVariantObj2, {'parameter','kRs', 'Value', estValues2.Estimate(2)});
 addcontent(optimVariantObj2, {'parameter','kRD1','Value', estValues2.Estimate(3)});
 addcontent(optimVariantObj2, {'parameter','kGd', 'Value', estValues2.Estimate(4)});

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

 optimVariantObj.Active = false;
 optimVariantObj2.Active = true;
 sdEst2 = sbiosimulate(modelObj);

Сравните с предыдущими результатами.

 [t2, GaFracEst2] = selectbyname(sdEst2, 'GaFrac');
 sdEst2Resampled = resample(sdEst2, tExpt, 'pchip');
 [~, GaFracEst2Resampled] = selectbyname(sdEst2Resampled, 'GaFrac');
 sse2 = norm(GaFracExpt - GaFracEst2Resampled)^2;
 rSquare2 = 1-sse2/sst;
 figure(fh);
 plot(t2, GaFracEst2, 'g-');
 legendText{end+1} = sprintf('4 Constants Changed, R^2 = %4.2f', rSquare2);
 legend(legendText{:});

С оценкой параметра, свободной варьироваться четыре параметра, подгонка к экспериментальным данным улучшилась далее. Отображенные итерации оптимизации показывают, что целевая функция уменьшилась, и значение R-квадрата увеличилось.

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

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

Заключение

Этот пример ввел различные аспекты функциональности SimBiology для построения моделей, симуляции и анализа с помощью модели сигнализации белка G.