Ссылка

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

Цели

Фон

В дрожжах Saccharomyces cerevisiae передача сигналов G-белка в реакции спаривания является хорошо охарактеризованным путем передачи сигнала. Феромон, секретируемый альфа-клетками, активирует рецептор альфа-фактора, связанный с G-белком (Ste2p), в клетках «а», что приводит к разнообразным клеточным реакциям, включая остановку клеточного цикла и синтез новых белков. 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 с названием «Гетеротримерный G-белок 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);

Назначьте ParameterVariiveNames в объекте кинетического закона. Это сопоставляет ParameterVariiveNames в объекте кинетического закона, чтобы можно было определить скорость реакции.

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)

Скорость реакции для первой реакции настроена.

 reactionObj1.ReactionRate

ans =

    'kRL*L*R - kRLm*RL'

Завершение создания модели Wild-Type

Чтобы создать модель штамма дикого типа (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';

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

 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 на 600. Кроме того, не регистрировать данные для лиганда 'L' (modelObj.Species (1)), потому что он принимает значения, которые на порядки выше, чем другие виды. Это делает визуализацию вида на участке более удобной. Для этого определите параметры «» Состояние «» для включения всех видов, кроме «» 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). Эта скорость определяется значением параметра скорости kGd. Мутантный штамм можно представить с помощью объекта Variant. Вариант SimBiology хранит альтернативные значения для одного или нескольких свойств модели SimBiology, такие как InitityAmount вида или Value параметра.

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

 variantObj = addvariant(modelObj, 'mutant');

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

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

Моделирование мутантного пути и построение графика результатов

Смоделировать модель с помощью объекта мутантного варианта. В процессе моделирования к параметру kGd применяется значение 0,004. Возврат результатов моделирования в объект 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-белка значительно ниже в мутантном штамме по сравнению с диким типом (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. Каждый объект содержит данные одного прогона сканирования параметров.

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

 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)

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

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

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

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

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

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

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

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

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

Оценка параметров игнорирует опцию анализа чувствительности, если она включена в конфигурационном наборе. Во избежание предупреждений отключите функцию «Анализ ошибок» в параметрах решателя.

 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-белка.