酵母异质三聚体G蛋白循环的参数扫描、参数估计和敏感性分析
这个例子展示了如何使用来自文献的途径在SimBiology®中构建、模拟和分析一个模型。
参考
酵母异质三聚体G蛋白循环的定量表征。伊头武，北野博明，梅尔文·西蒙。PNAS(2003)第100卷，10764-10769。
目标
为酵母TMY101(wt)菌株创建一个模型，显示g蛋白失活的野生型(催化)速率。
为TMY111(mutt)菌株创建一个变种，显示g蛋白失活的突变(非催化)率。
模拟并存储来自两个模型的数据。
比较野生型途径、突变型途径和实验数据之间g蛋白激活的时间进程。
执行参数扫描，以确定更改参数值对感兴趣的物种的影响。
利用实验数据估计模型参数值。
执行敏感性分析以确定模型参数对感兴趣的物种有哪些影响，以及影响到什么程度。
背景
在酵母中，G蛋白在交配反应中的信号转导是一个特征明确的信号转导途径。细胞分泌的信息素激活“a”细胞中的g蛋白偶联因子受体(Ste2p)，从而导致各种细胞反应，包括细胞周期阻滞和新蛋白质的合成。G蛋白和G蛋白偶联受体(GPCRs)是制药行业药物发现的重点。许多已上市的药物以GPCRs为靶点——一些例子包括减少胃酸(雷尼替丁，以组胺H2受体为靶点)、偏头痛(舒马曲坦，以血清素受体亚型为靶点)、精神分裂症(奥氮平，以血清素和多巴胺受体为靶点)和过敏(地氯雷他定，以组胺受体为靶点)。此外，一些估计表明，GPCRs是40%的药物发现努力的目标重点。一种方法是建立GPCR信号通路模型，分析和预测下游效应和相关通路的效应。这个例子检验了酵母信息素反应途径中G蛋白循环的模型构建、模拟和分析。
这张图是用来模拟酵母G蛋白循环的概念框架的图形表示。
使用下列缩写:
L =配体
R =受体
Gd = G - GDP
Gbg =游离的G - γ水平
Ga = G alpha - GTP
非活性异三聚体G蛋白(包含G-和G-)
Null =源或接收
Sst2为G蛋白调节因子(RGS) Sst2p
通路的反应
如图所示，这个循环可以浓缩成一组生化反应:
受体-配体相互作用(可逆反应)
+ r < - >
2)异质三聚体g蛋白形成
Gd + Gbg - >g
3) g -蛋白活化-下面注意，RL出现在方程的两边，因为RL是反应的修饰剂或催化剂。RL既不产生也不消耗。
RL + G -> Ga + Gbg + RL
4)受体合成与降解(以可逆反应表示降解与合成)
R < - > null
5) Receptor-ligand退化
RL - >零
6) g蛋白失活，在野生型菌株TMY101中由Sst2p催化，在突变型菌株TMY111中由SST2基因破坏未催化。
Ga - > Gd
所有的值已经转换为分子的种类数量，分子/秒，或1/秒的速率参数。
建立野生型通路的SimBiology®模型
创建一个名为“异质三聚体G蛋白wt”的SimBiology模型对象。
modelObj = sbiomodel ('异质三聚体G蛋白wt'）;
加入受体-配体相互作用(可逆反应)。
reactionObj1 = addreaction (modelObj,' l + r <-> rl '，．..“名字”，“Receptor-ligand互动”）;
对这个反应使用“质量作用”动力学定律。该模型是利用所有反应的质量作用动力学建立的。
kineticlawObj1 = addkineticlaw (reactionObj1,“MassAction”）;
添加正向和反向速率参数。
addparameter (modelObj“kRL”3.32 e-18);addparameter (modelObj“kRLm”, 0.01);
在动力学定律对象中分配parametervariablename。这将参数变量映射到动力学定律对象中的参数变量名，从而可以确定反应速率。
kineticlawObj1。ParameterVariableNames = {“kRL”，“kRLm”};
SimBiology自动为每个参与反应的物种创建物种对象。设定这些物种的初始数量。
%为“L”设置初始金额modelObj.Reactions (1) .Reactants(1)。InitialAmount = 6.022 e17;为“R”设置初始数量.Reactants modelObj.Reactions(1)(2)。InitialAmount = 10000.0;%保留“RL”的初始值为默认值(0.0)
第一个反应的反应率现在已经配置好了。
reactionObj1。ReactionRate
ans = 'kRL*L*R - kRLm*RL'
完成Wild-Type模型
创建野生型应变(TMY101)模型，添加其余的反应和参数，为每个反应创建动力学定律对象，并为动力学定律指定参数变量。
加入并配置异质三聚体g蛋白形成的反应。
reactionObj2 = addreaction (modelObj,'Gd + Gbg -> G'，．..“名字”，“G蛋白复合物形成”）;kineticlawObj2 = addkineticlaw (reactionObj2,“MassAction”）;addparameter (modelObj“kG1”, 1.0);kineticlawObj2。ParameterVariableNames =“kG1”；%为“Gd”设置初始数量modelObj.Reactions (2) .Reactants(1)。InitialAmount = 3000;%为“Gbg”设置初始数量modelObj.Reactions (2) .Reactants(2)。InitialAmount = 3000;%为“G”设置初始金额modelObj.Reactions (2)下载188bet金宝搏 . product(1)。InitialAmount = 7000;
添加并配置g蛋白活化反应。
reactionObj3 = addreaction (modelObj,'G + RL -> Ga + Gbg + RL'，．..“名字”，“G蛋白激活”）;kineticlawObj3 = addkineticlaw (reactionObj3,“MassAction”）;addparameter (modelObj“kGa”1.0 e-5);kineticlawObj3。ParameterVariableNames =“kGa”；为“Ga”设置初始数量modelObj.Reactions (3)下载188bet金宝搏 . product(1)。InitialAmount = 0.0;
添加和配置反应，以进行受体的合成和降解。
reactionObj4 = addreaction (modelObj,“R < - >空”，．..“名字”，“合成R /退化”）;kineticlawObj4 = addkineticlaw (reactionObj4,“MassAction”）;addparameter (modelObj“kRdo”, 4.0的军医);addparameter (modelObj“kr”, 4.0);kineticlawObj4。ParameterVariableNames = {“kRdo”，“kr”};
添加和配置反应的受体-配体降解。
reactionObj5 = addreaction (modelObj,“RL - >空”，“名字”，“RL退化”）;kineticlawObj5 = addkineticlaw (reactionObj5,“MassAction”）;addparameter (modelObj“kRD1”, 0.0040);kineticlawObj5。ParameterVariableNames =“kRD1”；
添加并配置g蛋白失活反应。
reactionObj6 = addreaction (modelObj,“Ga - > Gd”，“名字”，“Gprotein失活的）;kineticlawObj6 = addkineticlaw (reactionObj6,“MassAction”）;addparameter (modelObj“kGd”, 0.11);kineticlawObj6。ParameterVariableNames =“kGd”；
检查所有反应的反应速率。
get (modelObj。反应,{“反应”，“ReactionRate”}）
ans = 6 x2单元阵列{' L + R < - > RL”}{“kRL * L * R - kRLm * RL”}{“Gd + Gbg - > G”}{kG1 * * Gbg Gd的}{' G + RL - > Ga + Gbg + RL '}{“kGa * G * RL”}{' R < - >空'}{“kRdo * R - kr”}{RL - >空的}{“kRD1 * RL”}{Ga - > Gd的}{kGd * Ga的}
模拟野生型模型并绘制结果
为了记录物种Ga的快速上升和随后的下降，模拟模型600s并存储结果。
将默认配置设置对象的“停止时间”从10s (simulationTime)修改为600s。此外，不要记录配体'L' (modelObj.Species(1))的数据，因为它的值比其他物种高几个数量级。这使得物种在图中可视化更加方便。为了实现这一点，定义StatesToLog包括除'L'以外的所有物种。
configsetObj = getconfigset (modelObj);configsetObj。StopTime = 600;configsetObj.SolverOptions.AbsoluteTolerance = 1. e-9;configsetObj.RuntimeOptions.StatesToLog =．..sbioselect (modelObj“类型”，“物种”，“在那里”，“名字”，“~ = '，“L”）;
模拟模型并将结果返回给三个变量'time'、'data'和'names'。
[time, data, names] = sbiosimulate(modelObj);
图数据。
情节(时间、数据);传奇(名称,“位置”，“NorthEastOutside”）;包含(的时间(秒)）;ylabel (物种数量的）;网格在；
为突变株创建一个模型变体
突变株的g蛋白循环模型不同于活性g蛋白(Ga)失活的速率。该速率由速率参数kGd的值决定。您可以使用Variant对象表示突变菌株。SimBiology Variant存储SimBiology模型的一个或多个属性的替代值，例如物种的InitialAmount或参数的Value。
向模型添加一个名为“mutant”的变体。
variantObj = addvariant (modelObj,“变异”）;
向变体中添加内容，以为参数kGd指定0.004的备用值。
addcontent (variantObj, {“参数”，“kGd”，“价值”, 0.004});
模拟突变通路并绘制结果
使用变异变量对象模拟模型。这将0.004的值应用于模拟期间的参数kGd。在SimData对象中返回模拟结果。除了存储SimBiology模拟数据外，SimData对象还提供了数据访问、绘图和分析的方法。
将变异变量对象的Active属性设置为true并模拟。
variantObj。积极= true;mutantData = sbiosimulate (modelObj);
使用虚线绘制数据。另请参阅sbioplot方便地绘制SimData对象。
情节(mutantData。时间,mutantData。数据,“线型”，“——”）;传奇(mutantData。DataNames,“位置”，“NorthEastOutside”）;包含(的时间(秒)）;ylabel (物种数量的）;网格在；
比较活性g蛋白物种(Ga)在野生型和突变途径中的行为。
GaIndex = strcmp(名称,“遗传算法”）;%索引野生型结果[tmut, xmut] = selectbyname(mutantData, xmut)“遗传算法”）;plot(time, data(:，GaIndex)， tmut, xmut，“——”）;包含(的时间(秒)）;ylabel (物种数量的）;传奇({“Ga (wt)”，“Ga(突变)”}，“位置”，“NorthEastOutside”）;网格在；
执行参数扫描
与野生型相比，突变株的g蛋白失活率要低得多(kGd = 0.004 vs . 0.11)，这解释了在上述比较中观察到的随着时间的推移g蛋白(Ga)活性较高的原因。要更详细地了解kGd的变化如何影响Ga的级别，请执行几个模拟的参数扫描，其中kGd的值在一个范围内变化。下面的示例演示了对5个kGd值的参数扫描;要增加迭代次数，请更改参数中的值linspace下面的函数。
生成5个均匀间隔的kGd值，范围从0.001到0.15。
kGdValues = linspace(1e- 3,0.15, 5);
将参数扫描的结果存储在SimData对象数组中。初始化一个变量以保存该数组。
scanData = [];
为加速仿真准备模型。
sbioaccelerate (modelObj);
循环kGdValues并对每个值执行模拟。使用模型上的突变变量来修改仿真过程中使用的kGd值。
为kGd = kGdValues%在变体中设置所需的kGd值。variantObj。内容{1}{4}= kGd;%模拟模型，将结果存储在SimData对象中。sd = sbiosimulate (modelObj);scanData = [scanData;sd);结束
scanData现在是一个包含五个元素的SimData对象数组。每个对象包含参数扫描中一次运行的数据。
从SimData对象数组中提取遗传算法的时间进程，并在单个轴上绘图。下面的代码一步一步地构造情节;另外,看到sbioplot和sbiosubplot．
[scan, xscan] = selectbyname(scanata, xscan)“遗传算法”）;跳频=图;持有在；为C = 1:5 plot(scan{C}， xscan{C});str = sprintf (' k = %5.3f', kGdValues (c));文本(tscan {c}(结束),xscan {c}(结束),str);结束%注释情节。轴(gca (fh),“广场”）;标题(“改变kGd的值:对Ga的影响”）;包含(的时间(秒)）;ylabel (物种数量的）;网格在；持有从；
参数估计-背景
在对生物系统建模时，常常需要包括其数值未知或仅粗略知道的参数。如果系统中有一个或多个物种的实验数据，则可以通过改变它们并寻找那些能使模型模拟结果与实验数据最匹配的值来估计这些参数的值。
在本节的示例中，我们将探讨在试图将G蛋白模型与实验数据拟合的背景下的参数估计功能。
参数估计-比较模型结果与实验数据
在实验数据方面，参考文献图5中包含了G蛋白活性组分的时间历程。
存储实验时间和状态数据。
text = [0 10 30 60 110 300 450 600]';GaFracExpt = [0 0.35 0.4 0.36 0.39 0.33 0.24 0.17 0.2]';数据= groupedData(table(textexpt, 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
更改配置集上的RuntimeOptions以记录GaFrac。
configsetObj.RuntimeOptions.StatesToLog = GaFracObj;
使突变体失效。
variantObj。积极= false;
模拟模型，将结果存储在SimData对象中。
sdWild = sbiosimulate (modelObj);
为“GaFrac”获取数据，以便稍后在绘图中使用。
[twid, GaFracWild] = selectbyname(swid, swid, swid, swid)“GaFrac”）;
将仿真结果重新采样到实验时间向量上。
sdWildResampled = resample(sdWild, textexpt，“pchip”）;
获得“Ga”物种的重新采样数据。
[~， GaFracWildResampled] = selectbyname(sdWildResampled，“GaFrac”）;
计算r平方值，测量模拟数据和实验数据之间的拟合。
sst = norm(GaFracExpt - mean(GaFracExpt))^2;sse = norm(GaFracExpt - GaFracWildResampled)^2;rSquare = 1-sse / sst;
将仿真结果与遗传算法的实验数据进行对比。
跳频=图;情节(tExpt GaFracExpt,“罗”）;legendText = {“实验”};标题(“适合GaFrac实验数据”）;包含(的时间(秒)）;ylabel (物种数量的）;持有在；情节(tWild GaFracWild);legendText{结束+ 1}= sprintf ('原始，R^2 = %4.2f', rSquare);传奇(legendText {:});网格在；
参数估计-估计单个模型参数
由参数扫描可知，参数kGd的取值对Ga物种的时间进程有显著影响。让我们看看是否可以通过改变kGd的值来改善模型结果与实验数据的拟合。
根据实验数据进行参数估计，优化kGd值。在优化过程中绘制关于迭代的信息，最多15个迭代。
paramToEst = estimatedInfo (“kGd”）;kGdObj = sbioselect (modelObj,“名字”，“kGd”）;选择= optimset (“PlotFcns”@optimplotfval,“麦克斯特”15);result1 = sbiofit(modelObj, data，“GaFrac = GaFracExpt”paramToEst,．..[],“fminsearch”、选择);estValues1 = result1。编写此表达式ParameterEstimates
估计StandardError estValues1 = 1 x3表名称 _______ ________ _____________ {' kGd} 0.12142 - 0.0018542
将kGd的估价值存储在一个新的模型变量中。
optimVariantObj = addvariant (modelObj,“优化kGd”）;addcontent (optimVariantObj, {“参数”，“kGd”，“价值”, estValues1.Estimate});
激活新的变异，使“突变”的变异失活。
optimVariantObj。积极= true;mutantVariantObj = getvariant (modelObj,“变异”）;mutantVariantObj。积极= false;
利用kGd的估计值模拟模型。
sdEst1 = sbiosimulate (modelObj);
绘制GaFrac的数据并与之前的结果进行比较。
[t1, GaFracEst1] = selectbyname(sdEst1，“GaFrac”）;sdEst1Resampled = resample(sdEst1, textexpt，“pchip”）;[~， GaFracEst1Resampled] = selectbyname(sdEst1Resampled，“GaFrac”）;sse1 = norm(GaFracExpt - GaFracEst1Resampled)^2;rSquare1 = 1-sse1 / sst;图(跳频);情节(t1, GaFracEst1,“m -”）;legendText{结束+ 1}= sprintf ('kGd Changed, R^2 = %4.2f', rSquare1);传奇(legendText {:});
从R-square值可以看出，新的kGd估价值与实验数据的拟合略好一些。如果kGd的原始值只是一个粗略的估计，我们可以将这些结果解释为对原始估计的确认或对其的改进。
敏感性分析-背景
到目前为止，我们一直对活性G蛋白Ga的动态行为感兴趣。参数扫描显示，该物种显著地受到控制g蛋白失活的速率常数kGd的影响。通过参数估计，我们发现通过优化kGd的值，我们能够更好地拟合遗传算法的实验时间过程。
一个很自然的问题是，模型的其他参数会影响Ga级别，这些影响的程度是多少?敏感性分析允许您通过计算一个或多个物种(“输出”)相对于模型参数值或物种初始条件(“输入因子”)的随时间变化的导数来回答这些问题。
灵敏度分析-计算灵敏度
计算遗传算法对模型中各参数的灵敏度。将灵敏度完全规范化，以便相互比较。
使模型上的突变变量对象失效，以便在其初始值处使用kGd计算灵敏度。
optimVariantObj。积极= false;
在模型configset中建立灵敏度计算。
在求解器选项中打开灵敏度分析。configsetObj.SolverOptions.SensitivityAnalysis = true;%配置灵敏度分析的灵敏度输出和输入。sensitivityOpt = configsetObj.SensitivityAnalysisOptions;GaObj = sbioselect (modelObj,“类型”，“物种”，“名字”，“遗传算法”）;sensitivityOpt。输出= GaObj;params = sbioselect (modelObj,“类型”，“参数”，“在那里”，“名字”，“~ = '，“GaFrac”）;sensitivityOpt。输入=参数;sensitivityOpt。归一化=“全部”；
模拟模型。
sdSens = sbiosimulate (modelObj);
从SimData对象中提取灵敏度数据并绘制计算出的灵敏度。
[t, R, sensOutputs, sensInputs] = getsensmatrix(sdSens);R =紧缩(R);图;情节(t, R);标题(“遗传算法对各种参数的归一化灵敏度”）;包含(的时间(秒)）;ylabel (“敏感”）;传奇(sensInputs“位置”，“NorthEastOutside”）;网格在；
参数估计-估计多个参数
结果表明，Ga不仅对参数kGd敏感，而且对kGa、kRs和kRD1也敏感。(在情节中，其他敏感因素几乎是零。)改变这四个参数可以使其更符合实验数据。
估计这四个参数以匹配目标数据。使用前面配置的优化选项和模型中的当前参数值作为优化的起点。
选择参数kGa、kRs、kRD1、kGd进行估计。
paramsToEst = estimatedInfo ({“kGa”，“kr”，“kRD1”，“kGd”}）；
如果在配置集中启用了敏感性分析选项，参数估计将忽略该选项。在求解器选项中关闭灵敏度分析以避免警告。
configsetObj.SolverOptions.SensitivityAnalysis = false;result2 = sbiofit(modelObj，数据，“GaFrac = GaFracExpt”paramsToEst,．..[],“fminsearch”、选择);estValues2 = result2。ParameterEstimates
估计StandardError estValues2 = 4 x3表名称 ________ __________ _____________ {' kGa} 9.0068 3.0162 e-06 e-06{“kr”}4.549 - 11.867{‘kRD1} 0.0031018 0.0027647{‘kGd} 0.12381 - 0.053593
将这四个参数的估价值存储在一个新的模型变量中。
optimVariantObj2 = addvariant (modelObj,“四个参数优化”）;addcontent (optimVariantObj2, {“参数”，“kGa”，“价值”estValues2.Estimate (1)});addcontent (optimVariantObj2, {“参数”，“kr”，“价值”estValues2.Estimate (2)});addcontent (optimVariantObj2, {“参数”，“kRD1”，“价值”estValues2.Estimate (3)});addcontent (optimVariantObj2, {“参数”，“kGd”，“价值”, estValues2.Estimate (4)});
现在用新估计的参数值模拟模型。
optimVariantObj。积极= false;optimVariantObj2。积极= true;sdEst2 = sbiosimulate (modelObj);
与之前的结果进行比较。
[t2, GaFracEst2] = sdEst2，“GaFrac”）;sdEst2Resampled = resample(sdEst2, textexpt，“pchip”）;[~， GaFracEst2Resampled] = selectbyname(sdEst2Resampled，“GaFrac”）;sse2 = norm(GaFracExpt - GaFracEst2Resampled)^2;rSquare2 = 1-sse2 / sst;图(跳频);情节(t2, GaFracEst2“g -”）;legendText{结束+ 1}= sprintf ('4常量改变，R^2 = %4.2f', rSquare2);传奇(legendText {:});
参数估计可任意改变四个参数，进一步提高了对实验数据的拟合性。所显示的优化迭代结果表明，目标函数减小，r平方值增大。
注意，这里执行的四参数估计可能与生物学相关，也可能不相关，仅用于说明目的。
还请注意，将估计结果存储在变体中可以很容易地在模拟模型的不同版本之间来回切换。此时有四个版本:原始版本、突变版本和基于参数估计结果的两个版本。
结论
这个例子介绍了SimBiology使用G蛋白信号模型构建、模拟和分析模型的各个方面功能。
你也可以从以下列表中选择一个网站:
