酵母异三聚体G蛋白循环的参数扫描、参数估计和敏感性分析
本例展示了如何使用文献中的路径在SimBiology®中构建、模拟和分析模型。
参考
酵母异三聚体G蛋白循环的定量表征。伊头牧，北野博明，梅尔文·i·西蒙。PNAS(2003)第100卷，10764-10769。
目标
为酵母TMY101(wt)菌株创建一个模型，显示野生型(催化)g蛋白失活率。
为TMY111(mut)菌株创建一个变体，显示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 α和G β γ)
Null =源或接收器
Sst2为G蛋白调节因子(RGS) Sst2p
通路的反应
如图所示，循环过程可以浓缩为一组生化反应:
1)受体-配体相互作用(可逆反应)
L + r < - > rl
2)异三聚体g蛋白的形成
Gd + Gbg -> G
3) g蛋白活化-请注意，RL出现在方程两边，因为RL是反应的修饰剂或催化剂。这个反应既不产生RL也不消耗RL。
RL + G -> Ga + Gbg + RL
4)受体的合成和降解(按可逆反应处理，代表降解和合成)
R < - > null
5)受体-配体降解
RL -> null
6) g蛋白失活，在野生型菌株TMY101中由Sst2p催化，在突变型菌株TMY111中未催化，SST2基因被破坏。
Ga -> Gd
所有数值都转换为分子表示物种数量，分子/秒表示速率参数，或1/秒表示速率参数。
构建野生型通路的SimBiology®模型
创建名为“Heterotrimeric G Protein wt”的SimBiology模型对象。
modelObj = sbiommodel (异三聚体G蛋白wt）;
加入受体-配体相互作用(可逆反应)。
反动obj1 =地址(modelObj，' l + r <-> rl '，.．.“名字”，“Receptor-ligand互动”）;
使用“MassAction”动力学定律进行反应。该模型是使用所有反应的质量作用动力学建立的。
kineticlawObj1 = addkineticlaw(反动obj1，“MassAction”）;
添加正向速率和反向速率参数。
addparameter (modelObj“kRL”3.32 e-18);addparameter (modelObj“kRLm”, 0.01);
在动力学定律对象中分配ParameterVariableNames。这将在动力学定律对象中将ParameterVariables映射到ParameterVariableNames，以便可以确定反应速率。
kineticlawObj1。ParameterVariableNames = {“kRL”，“kRLm”};
SimBiology自动为每个参与反应的物种创建物种对象。设置这些物种的初始数量。
为'L'设置初始金额modelObj.Reactions (1) .Reactants(1)。InitialAmount = 6.022E17;为'R'设置初始金额.Reactants modelObj.Reactions(1)(2)。InitialAmount = 10000.0;将'RL'的初始金额保留为默认值(0.0)
第一个反应的反应速率现在已经配置好了。
reactionObj1。ReactionRate
ans = 'kRL*L*R - kRLm*RL'
完成野生型模型
为了创建野生型应变(TMY101)模型，添加其余的反应和参数，为每个反应创建动力学定律对象，并为动力学定律分配参数变量。
添加并配置异三聚体g蛋白形成的反应。
反动obj2 =地址(modelObj，'Gd + Gbg -> G'，.．.“名字”，“G蛋白复合物的形成”）;kineticlawObj2 = addkineticlaw(反动obj2，“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蛋白激活反应。
反动obj3 =地址(modelObj，'G + RL -> Ga + Gbg + RL'，.．.“名字”，“G蛋白激活”）;kineticlawObj3 = addkineticlaw(反动obj3，“MassAction”）;addparameter (modelObj“kGa”1.0 e-5);kineticlawObj3。ParameterVariableNames =“kGa”；%为“Ga”设置初始金额modelObj.Reactions (3)下载188bet金宝搏 . product(1)。InitialAmount = 0.0;
添加和配置反应的受体合成和降解。
反动obj4 =地址(modelObj，'R <-> null'，.．.“名字”，“合成R /退化”）;kineticlawObj4 = addkineticlaw(反动obj4，“MassAction”）;addparameter (modelObj“kRdo”, 4.0的军医);addparameter (modelObj“kr”, 4.0);kineticlawObj4。ParameterVariableNames = {“kRdo”，“kr”};
添加并配置受体-配体降解反应。
反动obj5 =地址(modelObj，'RL -> null'，“名字”，“RL退化”）;kineticlawObj5 = addkineticlaw(反动obj5，“MassAction”）;addparameter (modelObj“kRD1”, 0.0040);kineticlawObj5。ParameterVariableNames =“kRD1”；
添加并配置g蛋白失活反应。
反动obj6 =地址(modelObj，“Ga -> Gd”，“名字”，“Gprotein失活的）;kineticlaobj6 = addkineticlaw(反动obj6，“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并存储结果。
将默认配置集对象的StopTime从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] = sbiosimulation (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。主动=真;mutantData = sbiosimulation (modelObj);
用虚线绘制数据。另请参阅sbioplot方便绘制SimData对象。
情节(mutantData。时间,mutantData。数据,“线型”，“——”）;传奇(mutantData。DataNames,“位置”，“NorthEastOutside”）;包含(的时间(秒)）;ylabel (物种数量的）;网格在；
比较活跃的g蛋白物种(Ga)在野生型和突变途径中的行为。
GaIndex = strcmp(名称，“遗传算法”）;%野生型结果索引[tmut, xmut] = selectbyname(mutantData，“遗传算法”）;plot(time, data(:，GaIndex)， tmut, xmut，“——”）;包含(的时间(秒)）;ylabel (物种数量的）;传奇({“Ga (wt)”，“Ga(突变)”}，“位置”，“NorthEastOutside”）;网格在；
执行参数扫描
与野生型相比，突变株的g蛋白失活率要低得多(kGd = 0.004 vs kGd = 0.11)，这解释了在上述比较中观察到的随着时间推移，g蛋白(Ga)的活性水平较高。为了更详细地了解kGd的变化如何影响Ga的水平，请执行几个模拟的参数扫描，其中kGd的值在一个值范围内变化。下面的示例演示了对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 = sbiosimulation (modelObj);丑闻=[丑闻;sd);结束
scanData现在是SimData对象的五个元素数组。每个对象包含参数扫描中一次运行的数据。
从SimData对象数组中提取Ga的时间路线，并在单轴上绘制。下面的代码逐步构建情节;另外,看到sbioplot和sbiosubplot．
[tscan, xscan] = selectbyname(scanata, xscan)“遗传算法”）;Fh =图;持有在；为C = 1:5 plot(scan{C}， xscan{C});STR = sprintf(' k = %5.3f', kGdValues (c));Text (tscan{c}(end)， xscan{c}(end)， str);结束注释情节。轴(gca (fh),“广场”）;标题(“改变kGd值:对Ga的影响”）;包含(的时间(秒)）;ylabel (物种数量的）;网格在；持有从；
参数估计-背景
在对生物系统建模时，通常需要包括数值未知或仅大致知道的参数。如果系统中有一个或多个物种的实验数据，可以通过改变这些参数的值来估计这些参数的值，并寻找那些导致模型模拟结果与实验数据之间最佳拟合的值。
在本节示例中，我们将探索在尝试将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的绝对量，而是使用一个非常数参数和一个repeatedAssignment规则将这个分数添加到模型中。
GaFracObj = modelObj.addparameter(“GaFrac”，“ConstantValue”, 0);GaFracRule = modelObj.addrule('GaFrac = Ga / (Ga + G + Gd)'，“repeatedAssignment”）
GaFracRule = SimBiology规则数组索引:RuleType: Rule: 1 repeatedAssignment GaFrac = Ga / (Ga + G + Gd)
将配置集上的RuntimeOptions更改为记录GaFrac日志。
configsetObj.RuntimeOptions.StatesToLog = GaFracObj;
使突变变体失效。
variantObj。主动= false;
模拟模型，将结果存储在SimData对象中。
sdWild = sbiosimulation (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 =图;情节(tExpt GaFracExpt,“罗”）;legendText = {“实验”};标题(“符合GaFrac的实验数据”）;包含(的时间(秒)）;ylabel (物种数量的）;持有在；情节(tWild GaFracWild);legendText{end+1} = sprintf('原始，R^2 = %4.2f', rSquare);传奇(legendText {:});网格在；
参数估计-估计单个模型参数
从参数扫描可以看出，参数kGd的取值对Ga物种的时间进程有显著影响。让我们看看是否可以通过改变kGd的值来改善模型结果与实验数据的拟合。
根据实验数据进行参数估计，优化kGd值。在优化过程中绘制关于迭代的信息，最多可绘制15个迭代。
paramToEst = estimatedInfo(“kGd”）;kGdObj = sbioselect(modelObj，“名字”，“kGd”）;Opt = optimset(“PlotFcns”@optimplotfval,“麦克斯特”15);result1 = sbiofit(modelObj, data，'GaFrac = GaFracExpt'paramToEst,.．.[],“fminsearch”、选择);estValues1 = result1。ParameterEstimates
估计StandardError estValues1 = 1 x3表名称 _______ ________ _____________ {' kGd} 0.12142 - 0.0018548
在一个新的模型变量中存储kGd的估计值。
optimVariantObj = addvariant(modelObj，“优化kGd”）;addcontent (optimVariantObj, {“参数”，“kGd”，“价值”, estValues1.Estimate});
激活新变体，灭活“突变”变体。
optimVariantObj。主动=真;mutantVariantObj = getvariant(modelObj，“变异”）;mutantVariantObj。主动= false;
用kGd的估计值对模型进行模拟。
sdEst1 = sbiosimulation (modelObj);
绘制GaFrac的数据，并与之前的结果进行比较。
[t1, GaFracEst1] = selectbyname(sdEst1，“GaFrac”）;sdEst1Resampled = resample(sdEst1, tExpt，“pchip”）;[~， GaFracEst1Resampled] = selectbyname(sdEst1Resampled, sdEst1Resampled, sdEst1Resampled)“GaFrac”）;sse1 = norm(GaFracExpt - GaFracEst1Resampled)^2;rSquare1 = 1-sse1/sst;图(跳频);情节(t1, GaFracEst1,“m -”）;legendText{end+1} = sprintf('kGd已更改，R^2 = %4.2f', rSquare1);传奇(legendText {:});
从r平方值，我们看到，拟合的实验数据略好，kGd的新估计值。如果kGd的原始值只是一个粗略的估计，我们可以将这些结果解释为对原始估计的确认或对原始估计的改进。
敏感性分析-背景
到目前为止，我们一直对活性G蛋白Ga种的动态行为感兴趣。参数扫描显示，该物种明显受控制g蛋白失活的速率常数kGd值的影响。通过参数估计，我们发现通过优化kGd的值，我们能够更好地拟合Ga的实验时间过程。
一个自然的问题是，模型的其他参数影响Ga水平，这些影响的大小是什么?灵敏度分析允许您通过计算一个或多个物种(“输出”)相对于模型参数值或物种初始条件(“输入因素”)的时间依赖导数来回答这些问题。
灵敏度分析-计算灵敏度
计算Ga对模型中各参数的灵敏度。完全标准化敏感性，以便它们可以相互比较。
在模型上禁用突变的变体对象，以便用kGd的原始值计算灵敏度。
optimVariantObj。主动= false;
在模型configset中设置灵敏度计算。
在解算器选项中打开灵敏度分析。configsetObj.SolverOptions.SensitivityAnalysis = true;为灵敏度分析配置灵敏度输出和输入。敏感性选择= configsetobj .敏感性分析选项;GaObj = sbioselect(modelObj，“类型”，“物种”，“名字”，“遗传算法”）;sensitivityOpt。输出= GaObj;params = sbioselect(modelObj，“类型”，“参数”，“在那里”，“名字”，“~ = '，“GaFrac”）;sensitivityOpt。输入=参数;sensitivityOpt。归一化=“全部”；
模拟模型。
sdSens = sbiosimulation (modelObj);
从SimData对象中提取灵敏度数据并绘制计算出的灵敏度。
[t, R, sensInputs] = getsensmatrix(sdSens);R =挤压(R);图;情节(t, R);标题(“Ga对不同参数的归一化灵敏度”）;包含(的时间(秒)）;ylabel (“敏感”）;传奇(sensInputs“位置”，“NorthEastOutside”）;网格在；
参数估计-估计多个参数
结果表明，Ga不仅对kGd参数敏感，而且对kGa、kRs和kRD1参数也敏感。(从图表上看，其他敏感度几乎为零。)改变这四个参数可以使模型更好地拟合实验数据。
估计这四个参数以匹配目标数据。使用先前配置的优化选项和模型中的当前参数值作为优化的起点。
选择参数kGa、kRs、kRD1、kGd进行估计。
paramsToEst = estimatedInfo({“kGa”，“kr”，“kRD1”，“kGd”})；
如果在configset中启用了敏感性分析选项，参数估计将忽略该选项。关闭解算器选项中的敏感性分析以避免警告。
configsetObj.SolverOptions.SensitivityAnalysis = false;result2 = sbiofit(modelObj, data，'GaFrac = GaFracExpt'paramsToEst,.．.[],“fminsearch”、选择);estValues2 = result2。ParameterEstimates
估计StandardError estValues2 = 4 x3表名称 ________ __________ _____________ {' kGa} 9.0068 3.0249 e-06 e-06{“kr”}4.549 - 11.786{‘kRD1} 0.0031018 0.0027417{‘kGd} 0.12381 - 0.053702
将四个参数的估计值存储在一个新的模型变量中。
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。主动=真;sdEst2 = sbiosimulation (modelObj);
与之前的结果进行比较。
[t2, GaFracEst2] = selectbyname(sdEst2，“GaFrac”）;sdEst2Resampled = resample(sdEst2, tExpt，“pchip”）;[~， GaFracEst2Resampled] = selectbyname(sdEst2Resampled，“GaFrac”）;sse2 = norm(GaFracExpt - GaFracEst2Resampled)^2;rSquare2 = 1-sse2/sst;图(跳频);情节(t2, GaFracEst2“g -”）;legendText{end+1} = sprintf('4个常数改变，R^2 = %4.2f', rSquare2);传奇(legendText {:});
参数估计可以改变四个参数，进一步提高了与实验数据的拟合性。所显示的优化迭代结果表明，目标函数减小，r平方值增大。
请注意，这里执行的四个参数估计可能与生物学相关，也可能不相关，仅用于说明目的。
还要注意，将估计结果存储在变量中可以很容易地在模拟模型的不同版本之间来回切换。目前有四个版本:原始版本，突变版本，以及基于参数估计结果的两个版本。
结论
本示例介绍了SimBiology功能的各个方面，用于使用G蛋白信号模型构建、模拟和分析。
您也可以从以下列表中选择一个网站: