Ссылка

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

Цели

Фон

В дрожжах Saccharomyces cerevisiae передача сигналов G-белка в спариваемом ответе является хорошо охарактеризованным путём передачи сигнала. Феромон, секретируемые альфа- камеры, активирует рецептор (Ste2p) альфа-фактора, связанный с G-белком, в 'a' камер что приводит к различным камерам реакциям, включая остановку клеточного цикла и синтез новых белков. G-белки и G-белковые связанные рецепторы (GPCR) являются особым вниманием усилий по открытию лекарств в фармацевтической промышленности. Многие продаваемые лекарства нацелены на GPCR - некоторые примеры включают препараты для восстановления желудочной кислоты (ранитидин, мишеней гистамин H2 рецептор), мигрени (суматриптан, мишенью для подтипа серотонинового рецептора), шизофрении (оланзапин, мишеневые серотониновые и дофаминовые рецепторы) и Кроме того, некоторые оценки предполагают, что GPCR являются целевым особым вниманием 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 с именем 'Heterotrimeric G Protein wt'.

 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 объекта конфигурации модели по умолчанию с 10s (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'.

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

Постройте график данных.

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

Создание варианта модели для мутантного штамма

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

Добавьте к модели вариант с именем 'mutant'.

 variantObj = addvariant(modelObj, 'mutant');

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

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

Симуляция мутантного пути и графическое изображение результатов

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

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

 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-белка (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-белка намного ниже в мутантном штамме по сравнению с диким типом (кГд = 0,004 против кГд = 0,11), что объясняет более высокие уровни активированного G-белка (Га) с течением времени, наблюдаемые в вышеописанном сравнении. Для более подробного просмотра того, как изменение kGd влияет на уровни Ga, выполните скан параметра из нескольких симуляций, в которых значение kGd варьируется в области значений значений. Следующий пример иллюстрирует сканирование параметра над пятью значениями кГд; чтобы увеличить количество итераций, измените значения в аргументах для linspace функция ниже.

Сгенерируйте пять равномерно расположенных значений кГд в диапазоне от 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-белка к экспериментальным данным.

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

Для экспериментальных данных Fig. 5 ссылки статьи содержит график для фракции активного 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, добавьте эту дробь к модели с помощью непостоянного параметра и правила reputedAssignment.

 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 оказывает значительный эффект на график вида Ga. посмотрим, сможем ли мы улучшить подгонку результатов модели к экспериментальным данным путем изменения значения 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  

Сохраните предполагаемое значение кГд в новой модели Variant.

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

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

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

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

 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-квадрат мы видим, что подгонка к экспериментальным данным немного лучше с новым, оценочным значением кГд. Если бы исходное значение для кГд было только приблизительной оценкой, мы могли бы интерпретировать эти результаты либо как подтверждение исходной оценки, либо как улучшение по сравнению с ней.

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

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

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

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

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

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

 optimVariantObj.Active = false;

Настройте расчет чувствительности в конфигурации модели.

% 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, но и кГа, кРс и kRD1. (Другие чувствительности неотличимы от нуля на графике.) Изменение этих четырех параметров может сделать подгонку к экспериментальным данным лучше еще.

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

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

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

Оценка параметра будет игнорировать опцию анализа чувствительности, если она включена в наборе настроек. Отключите 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-белка.