扫描参数、参数估计和灵敏度分析在酵母Heterotrimeric G蛋白循环
这个例子展示了如何构建、模拟和分析模型在SimBiology®使用途径从文学。
参考
酵母的定量表征heterotrimeric G蛋白循环。总裁中西宏明Tau-Mu咦,北野,梅尔文。西蒙。PNAS(2003) 100卷,10764 - 10769。
的目标是
创建一个模型的酵母TMY101 (wt)显示了野生型菌株(催化)蛋白失活率。
创建变体TMY111(傻瓜)应变显示突变体(uncatalyzed)蛋白失活率。
从这两个模型模拟和存储数据。
野生型之间的比较timecourse蛋白激活通路,突变体途径和实验数据。
执行参数扫描来确定不同的价值的影响参数对感兴趣的一个物种。
使用实验数据估计模型参数值。
进行灵敏度分析,以确定哪些和模型参数在多大程度上影响一个物种的利益。
背景
酵母酿酒酵母,在交配的响应是一个G蛋白信号特征信号转导通路。信息素分泌α细胞激活g蛋白耦合alpha-factor受体(Ste2p) a细胞导致多种细胞反应包括细胞循环逮捕和合成新的蛋白质。G蛋白和G蛋白耦合受体(GPCRs)是药物发现工作的重点在制药行业。许多销售药物目标GPCRs——一些例子包括减少胃酸(雷尼替丁,目标组胺H2受体),偏头痛(舒马曲坦,目标5 -羟色胺受体亚型),精神分裂症(奥氮平,目标5 -羟色胺和多巴胺受体)和过敏(地氯雷他定,目标组胺受体)。此外,一些估计显示GPCRs目标40%的药物发现努力的焦点。一种方法是模型GPCR信号通路来分析和预测相关通路的下游效应和影响。这个示例检查模型建立、仿真和分析的G蛋白周期酵母信息素反应通路。
这个数字是一个概念性的框架用于模型的图形化表示酵母G蛋白循环。
以下使用缩写:
L =配体
R =受体
Gd = Gα- GDP
Gbgβγ=自由水平的G
Ga = Gα-三磷酸鸟苷
G =不活跃heterotrimeric蛋白(包含G Gα和βγ)
null =源或下沉
Sst2表示G蛋白调节器Sst2p(该公司)
通路的反应
如图,周期可以浓缩为一系列生化反应:
1)Receptor-Ligand交互(可逆反应)
L + R < - > RL
2)Heterotrimeric蛋白形成
Gd + Gbg - > G
3)蛋白激活-注意下面RL出现在方程两边因为RL修饰符或反应的催化剂。RL既不生产,也不被这个反应。
RL + G - > Ga + Gbg + RL
4)受体合成和降解(视为可逆反应代表降解和合成)
R < - > null
5)Receptor-ligand退化
RL - >零
6)蛋白失活,催化Sst2p野生型菌株TMY101和uncatalyzed突变株与中断TMY111 SST2基因。
Ga - > Gd
所有值已经转化为分子物种数量,和分子/秒或1 /秒速度参数。
构建一个SimBiology®野生型路径模型
创建一个SimBiology模型对象的名称“Heterotrimeric克蛋白质wt”。
modelObj = sbiomodel (“Heterotrimeric克蛋白质wt”);
添加Receptor-Ligand交互(可逆反应)。
reactionObj1 = addreaction (modelObj,“L + R < - > RL”,…“名字”,“Receptor-ligand互动”);
使用“MassAction”动力法的反应。该模型使用质量作用动力学反应。
kineticlawObj1 = addkineticlaw (reactionObj1,“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.022 e17;%设置初始金额' R '.Reactants modelObj.Reactions (1) (2)。InitialAmount = 10000.0;%离开初始金额“RL”默认值(0.0)
第一反应的ReactionRate现在已经配置。
reactionObj1.ReactionRate
ans = ' kRL * L * R - kRLm * RL”
完成了野生型模型
创建模型的野生型菌株(TMY101),添加其余的反应和参数,创建动态法律对象的反应,并分配参数变量的动力学规律。
添加和配置Heterotrimeric 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;
添加和配置蛋白活化的反应。
reactionObj3 = addreaction (modelObj,“G + RL - > Ga + Gbg + RL”,…“名字”,“G蛋白激活”);kineticlawObj3 = addkineticlaw (reactionObj3,“MassAction”);addparameter (modelObj“kGa”1.0 e-5);kineticlawObj3。ParameterVariableNames =“kGa”;%为“遗传算法”设定初始量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”};
添加和配置receptor-ligand降解的反应。
reactionObj5 = addreaction (modelObj,“RL - >空”,“名字”,“RL退化”);kineticlawObj5 = addkineticlaw (reactionObj5,“MassAction”);addparameter (modelObj“kRD1”,0.0040);kineticlawObj5。ParameterVariableNames =“kRD1”;
添加和配置反应蛋白失活。
reactionObj6 = addreaction (modelObj,“Ga - > Gd”,“名字”,“Gprotein失活的);kineticlawObj6 = addkineticlaw (reactionObj6,“MassAction”);addparameter (modelObj“kGd”,0.11);kineticlawObj6。ParameterVariableNames =“kGd”;
检查所有的ReactionRate反应。
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的}
模拟野生型模型并绘制结果
要注意的快速崛起和随后的衰落物种遗传算法,模拟模型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“类型”,“物种”,“在那里”,“名字”,“~ = ',“L”);
模拟模型,并将结果返回给三个变量“时间”,“数据”和“名字”。
[时间、数据名称]= sbiosimulate (modelObj);
图数据。
情节(时间、数据);传奇(名称,“位置”,“NorthEastOutside”);包含(的时间(秒));ylabel (物种数量的);网格在;
创建模型变异的突变株
突变菌株的蛋白周期模型不同的失活速率活性蛋白(Ga)。这个速度是由kGd率参数的值。你可以使用变体对象代表突变株。SimBiology变体商店替代值SimBiology模型的一个或多个属性,如InitialAmount一个物种或一个参数的值。
添加一个名为“突变”模型的变体。
variantObj = addvariant (modelObj,“变异”);
将内容添加到指定另一种变体kGd参数值为0.004。
addcontent (variantObj, {“参数”,“kGd”,“价值”,0.004});
模拟突变通路并绘制结果
模拟对象模型使用突变变体。这0.004的值适用于参数kGd在模拟。仿真结果在SimData对象返回。除了存储SimBiology仿真数据,SimData对象为数据访问提供方法,绘制和分析。
突变体变体对象的Active属性设置为真实和模拟。
variantObj。积极= true;mutantData = sbiosimulate (modelObj);
用虚线画出数据。另请参阅sbioplot为方便SimData对象的绘制。
情节(mutantData。时间、mutantData.Data“线型”,“——”);传奇(mutantData.DataNames“位置”,“NorthEastOutside”);包含(的时间(秒));ylabel (物种数量的);网格在;
比较活跃的行为g蛋白种类(Ga)野生型和突变体通路。
GaIndex = strcmp(名称,“遗传算法”);%指数野生型的结果[tmut, xmut] = selectbyname (mutantData,“遗传算法”);情节(时间、数据(:,GaIndex) tmut, xmut,“——”);包含(的时间(秒));ylabel (物种数量的);传奇({“Ga (wt)”,“Ga(突变)”},“位置”,“NorthEastOutside”);网格在;
执行参数扫描
蛋白失活率低得多的突变株与野生型(kGd = 0.004 vs kGd = 0.11),这就解释了更高水平的激活蛋白(Ga)随着时间的推移,观察到在上面的比较。进行更详细的观察kGd的变化如何影响水平的Ga,执行参数扫描的几个模拟kGd多样的价值在一个范围的值。下面的例子说明了一个参数扫描在五kGd值;增加迭代的数量,改变参数的值linspace下面的函数。
生成5个等间隔kGd值从0.001到0.15不等。
kGdValues = linspace (1 e - 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对象的数组。每个对象包含的数据从一个运行参数扫描。
提取timecourses SimData Ga的对象数组和情节在一个轴上。下面的代码构造情节循序渐进;另外,看到sbioplot和sbiosubplot。
[tscan, xscan] = selectbyname (scanData,“遗传算法”);跳频=图;持有在;为c = 1:5情节(tscan c {}, xscan {c});str = sprintf (“k = % 5.3 f ',kGdValues (c));文本(tscan {c}(结束),xscan {c}(结束),str);结束%注释情节。轴(gca (fh),“广场”);标题(的不同的价值kGd:影响遗传算法的);包含(的时间(秒));ylabel (物种数量的);网格在;持有从;
参数估计- - -背景
生物系统建模时,常常需要包括参数的数值是未知的或只大概知道。如果实验数据可用于一个或多个物种的系统,估计这些参数的值可以通过改变他们,寻找这些价值观导致之间的最佳匹配模型的模拟结果和实验数据。
在本节的例子中,我们探索的上下文中参数估计功能试图适应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(表(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:规则:1 repeatedAssignment GaFrac = Ga / (Ga + G + Gd)
更改配置设置日志GaFrac RuntimeOptions。
configsetObj.RuntimeOptions。StatesToLog = GaFracObj;
停用突变变体。
variantObj。积极= false;
模拟模型,将结果存储在一个SimData对象。
sdWild = sbiosimulate (modelObj);
得到的数据“GaFrac”在一块使用。
[tWild, GaFracWild] = selectbyname (sdWild,“GaFrac”);
重新取样仿真结果在实验时间向量。
sdWildResampled =重新取样(sdWild tExpt,“pchip”);
得到重新取样数据的物种“遗传算法”。
[~,GaFracWildResampled] = selectbyname (sdWildResampled,“GaFrac”);
计算r平方值测量模拟和实验数据之间的配合。
海温=规范(GaFracExpt -意味着(GaFracExpt)) ^ 2;上交所=规范(GaFracExpt GaFracWildResampled) ^ 2;rSquare = 1-sse / sst;
情节为Ga仿真结果与实验数据。
跳频=图;情节(tExpt GaFracExpt,“罗”);legendText = {“实验”};标题(“GaFrac适合实验数据”);包含(的时间(秒));ylabel (物种数量的);持有在;情节(tWild GaFracWild);legendText{结束+ 1}= sprintf (原来,R ^ 2 = % 4.2 f ',rSquare);传奇(legendText {:});网格在;
参数估计,估计模型参数
从参数扫描,我们看到,参数的值kGd显著影响物种的timecourse Ga。看看我们可以提高模型的适应结果的实验数据不同kGd的价值。
实验数据进行参数估计,优化kGd的价值。情节信息迭代优化过程,最多15个迭代。
paramToEst = estimatedInfo (“kGd”);kGdObj = sbioselect (modelObj,“名字”,“kGd”);选择= optimset (“PlotFcns”@optimplotfval,“麦克斯特”15);result1 = sbiofit编写此表达式(modelObj、数据“GaFrac = GaFracExpt”paramToEst,…[],“fminsearch”、选择);estValues1 = result1.ParameterEstimates
estValues1 = 1 x3表名字估计StandardError _________ ________ _________________ {‘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 =重新取样(sdEst1 tExpt,“pchip”);[~,GaFracEst1Resampled] = selectbyname (sdEst1Resampled,“GaFrac”);sse1 =规范(GaFracExpt GaFracEst1Resampled) ^ 2;rSquare1 = 1-sse1 / sst;图(跳频);情节(t1, GaFracEst1,“m -”);legendText{结束+ 1}= sprintf (“kGd改变,R ^ 2 = % 4.2 f ',rSquare1);传奇(legendText {:});
从r平方值,我们看到,适合实验数据与新略好,估计kGd的价值。如果原始值kGd只是一个粗略的估计,我们可以解释这些结果作为确认最初的估计或有所改善。
敏感性分析,背景
到目前为止,我们一直感兴趣的动态行为活跃的G蛋白,物种遗传算法。参数扫描显示,这个物种是影响速率常数的值kGd执政蛋白失活。使用参数估计,我们发现,通过优化kGd的价值,我们能够更好地适应一个实验性timecourse Ga。
一个自然的问题是,其他参数模型的影响遗传算法的水平,这些影响的大小是什么?灵敏度分析可以回答这些问题通过计算一个或多个物种的时间衍生品(“输出”)相对于模型参数值,或物种的初始条件(“输入因素”)。
敏感性分析,计算敏感性
计算的敏感性Ga对模型中的各种参数。规范化的敏感性完全可以相互比较。
停用突变变体对象与kGd模型以便敏感性计算在其原始值。
optimVariantObj。积极= false;
configset设置灵敏度计算的模型。
%打开SensitivityAnalysis解算器的选项。configsetObj.SolverOptions。SensitivityAnalysis = true;%配置灵敏度输出和输入的灵敏度分析。sensitivityOpt = configsetObj.SensitivityAnalysisOptions;GaObj = sbioselect (modelObj,“类型”,“物种”,“名字”,“遗传算法”);sensitivityOpt。输出= GaObj;params = sbioselect (modelObj,“类型”,“参数”,“在那里”,“名字”,“~ = ',“GaFrac”);sensitivityOpt。输入=参数;sensitivityOpt。归一化=“全部”;
模拟模型。
sdSens = sbiosimulate (modelObj);
从SimData提取敏感数据对象和绘制计算敏感性。
R (t)、sensOutputs sensInputs] = getsensmatrix (sdSens);R =紧缩(R);图;情节(t, R);标题(“规范化Ga对各种参数的敏感性);包含(的时间(秒));ylabel (“敏感”);传奇(sensInputs“位置”,“NorthEastOutside”);网格在;
参数估计,估计多个参数
这些结果表明,Ga kGd不仅是敏感参数,而且kGa,基米-雷克南,kRD1。(其他敏感的零阴谋。)不同这四个参数可以更好的适合的实验数据。
估计这四个参数,以匹配目标数据。使用先前配置的优化选择和当前模型中参数值作为优化的起点。
选择参数kGa、kr kRD1, kGd估计。
paramsToEst = estimatedInfo ({“kGa”,“kr”,“kRD1”,“kGd”});
参数估计会忽略了灵敏度分析选项如果configset中启用它。在解算器选项关掉SensitivityAnalysis避免警告。
configsetObj.SolverOptions。SensitivityAnalysis = false;result2 = sbiofit (modelObj、数据“GaFrac = GaFracExpt”paramsToEst,…[],“fminsearch”、选择);estValues2 = result2.ParameterEstimates
estValues2 = 4 x3表名估计StandardError ________ __________ _________________ {‘kGa} 9.0068 e-06 3.0162 e-06 4.549 - 11.867 {“kr”} {‘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] = selectbyname (sdEst2“GaFrac”);sdEst2Resampled =重新取样(sdEst2 tExpt,“pchip”);[~,GaFracEst2Resampled] = selectbyname (sdEst2Resampled,“GaFrac”);sse2 =规范(GaFracExpt GaFracEst2Resampled) ^ 2;rSquare2 = 1-sse2 / sst;图(跳频);情节(t2, GaFracEst2“g -”);legendText{结束+ 1}= sprintf (“4常数改变,R ^ 2 = % 4.2 f ',rSquare2);传奇(legendText {:});
参数估计自由改变四个参数,适合实验数据进一步改善了。显示的优化迭代显示目标函数下降和r平方值也增加了。
注意,这里执行的四个参数估计可能是也可能不是生物学上相关,仅为演示之用。
还请注意,将评估结果存储在变量便于之间来回切换模拟模型的不同版本。在这一点上有四个版本:最初,突变体,两个版本基于参数估计的结果。
结论
这个例子介绍的各个方面SimBiology构建功能模型,模拟和分析使用G蛋白信号的典范。
