酵母异源三聚体G蛋白周期中的参数扫描、参数估计和敏感性分析

开放脚本

这个例子展示了如何使用来自文献的途径在SimBiology®中构建、模拟和分析一个模型。

参考

酵母异络蛋白循环的定量表征。Tau-mu yi，Hiroaki Kitano，和Melvin I. Simon。PNAS（2003）卷。100,10764-10769。

目的

为酵母TMY101(wt)菌株创建一个模型，显示g蛋白失活的野生型(催化)速率。

为TMY111（mut）菌株创建一个显示G蛋白失活突变（未催化）率的变体。

模拟并存储来自两个模型的数据。

比较野生型途径、突变途径和实验数据之间G蛋白激活的时间过程。

执行参数扫描以确定在感兴趣的物种上改变参数的值的效果。

使用实验数据估计模型参数值。

进行敏感性分析，以确定模型参数对感兴趣物种的影响程度。

出身背景

在酿酒酵母中，交配反应中的G蛋白信号传导是一个具有良好特征的信号转导途径。α细胞分泌的信息素激活“a”细胞中的G蛋白偶联α因子受体（Ste2p），导致多种细胞反应，包括细胞周期阻滞和新蛋白质的合成。G蛋白和G蛋白偶联受体（GPCR）是制药行业药物研发工作的重点。许多上市药物以GPCR为目标-一些例子包括用于减少胃酸（雷尼替丁，靶向组胺H2受体）、偏头痛（舒马曲坦，靶向血清素受体亚型）、精神分裂症（奥氮平，靶向血清素和多巴胺受体）和过敏（地氯雷他定，靶向组胺受体）的药物。此外，一些估计表明，GPCR是40%药物发现工作的目标焦点。一种方法是对GPCR信号通路进行建模，以分析和预测下游效应和相关通路中的效应。本例检验酵母信息素反应途径中G蛋白循环的模型构建、模拟和分析。

该图是用于模拟酵母G蛋白质循环的概念框架的图形表示。

使用以下缩写：

L=配体
R=受体
gd = g alpha-gdp
Gbg=Gβγ的自由水平
ga = g alpha - gtp
G =无活性异质型G-蛋白（含有Gα和Gβ-Gamma）
null=源或接收器
Sst2表示G蛋白调节因子（RGS）Sst2p

通路的反应

如图所示，这个循环可以浓缩成一组生化反应:

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蛋白失活，在野生型菌株TMY101中由Sst2p催化，在突变型菌株TMY111中由SST2基因破坏未催化。

Ga->Gd

对于物种数量，所有值均已转换为分子；对于速率参数，所有值均已转换为分子/秒或1/秒。

建立野生型通路的SimBiology®模型

创建一个名为“异质三聚体G蛋白wt”的SimBiology模型对象。

modelObj = sbiomodel ('杂漆g蛋白质重量）;

添加受体-配体相互作用（可逆反应）。

reactionObj1 = addreaction (modelObj,“L+R<->RL”,...“姓名”,“受体-配体相互作用”）;

使用“批量”动力学法进行反应。该模型是使用质量动作动力学构建的所有反应。

kineticlawObj1 = addkineticlaw (reactionObj1,“按摩”）;

添加正向和反向速率参数。

addparameter（modelobj，“kRL”，3.32e-18）；添加参数（modelObj，“kRLm”, 0.01);

在动力学定律对象中指定ParameterVariableNames。这将ParameterVariables映射到动力学定律对象中的ParameterVariableNames，以便可以确定反应速率。

kineticlawObj1.ParameterVariableNames={“kRL”,“kRLm”}；

SimBiology自动为反应中的每个参与物种创建物种对象。设置这些物种的初始数量。

%为“L”设置初始金额模型对象反应（1）.反应物（1）.初始量=6.022E17；%设置“R”的初始金额模型对象反应（1）.反应物（2）.初始量=10000.0；%将“RL”的初始金额保留为默认值（0.0）

第一个反应的反应率现在已经配置好了。

反应速率

ans ='krl * l * r  -  krlm * rl'

完成Wild-Type模型

要创建野生型应变的模型（TMY101），请添加其余的反应和参数，为每个反应创建动力学定律对象，并为动力学法分配参数变量。

添加并配置异三聚体G蛋白形成反应。

reactionObj2=添加反应（modelObj，“Gd+Gbg->G”,...“姓名”,'G蛋白复杂的形成'）;kineticlawObj2 = addkineticlaw (reactionObj2,“按摩”); addparameter（modelObj，“kG1”，1.0）；kineticlawObj2.ParameterVariableNames=“kG1”；%为“Gd”设置初始数量modelobj.reactions（2）。重演（1）.initialamount = 3000;%为“Gbg”设置初始数量模型对象反应（2）.反应物（2）.初始量=3000；%设置“G”的初始金额模型对象反应（2）.产物（1）.初始量=70s manbetx 84500；

添加并配置g蛋白活化反应。

ressobj3 = addReacection（Modelobj，'G + RL -> Ga + Gbg + RL',...“姓名”,“G蛋白活化”）;KineticLawobj3 = addkineticlaw（ressionobj3，“按摩”); addparameter（modelObj，“克格勃”，1.0E-5）；运动学定律OBJ3.参数变量名称=“克格勃”；为“Ga”设置初始数量模型对象反应（3）.产物（1）.初始量=0.s manbetx 8450；

添加并配置受体合成和降解反应。

reactionObj4=添加反应（modelObj，“R<->null”,...“姓名”,'R合成/降解')；kineticlawObj4=添加kineticlaw（反应bj4，“按摩”); addparameter（modelObj，“克尔多”，4.0E-4）；添加参数（modelObj，“kRs”，4.0）；kineticlawObj4.ParameterVariableNames={“克尔多”,“kRs”}；

添加和配置反应的受体-配体降解。

reactionObj5=添加反应（modelObj，'RL->null',“姓名”,'rl劣化'）;KineticLawobj5 = Addkineticlaw（ressionobj5，“按摩”); addparameter（modelObj，“kRD1”, 0.0040);kineticlawObj5。ParameterVariableNames =“kRD1”；

添加并配置g蛋白失活反应。

ressobj6 = addreacrection（modelobj，“Ga - > Gd”,“姓名”,“G蛋白失活”）;KineticLawobj6 = Addkineticlaw（Astryovj6，“按摩”); addparameter（modelObj，'kgd'，0.11）;KineticLawobj6.ParametervariaBlenames ='kgd'；

检查所有反应的反应速率。

得到（modelobj.reaciptions，{'反应',“反应速率”}）

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的快速上升和随后的下降，模拟模型600秒并存储结果。

将默认配置集对象的停止时间从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，“类型”,“物种”,“哪里”,“姓名”,“~ = ',“我）;

模拟模型并将结果返回到三个变量“时间”、“数据”和“名称”。

[时间、数据、名称]=sbiosimulate（modelObj）；

绘制数据。

图（时间，数据）;传奇（名字，'地点',“东北朝”)；xlabel('时间（秒）')；伊莱贝尔(“物种数量”）;网格在；

为突变菌株创建模型变体

突变菌株的G蛋白循环模型在活性G蛋白（Ga）失活的速率上有所不同发生。此速率由速率参数kGd的值控制。您可以使用Variant对象表示突变菌株。SimBiology变体存储SimBiology模型的一个或多个属性的备选值，例如物种的初始数量或参数值。

向模型中添加名为“突变体”的变体。

variantObj=添加变量（modelObj，'突变'）;

向变体中添加内容，以为参数kGd指定0.004的备用值。

addcontent（variantObj{“参数”,'kgd',“价值”，0.004}）;

模拟突变途径并绘制结果

使用突变体变体对象模拟模型。这在仿真期间将0.004的值应用于参数KGD。返回Simdata对象中的模拟结果。除了存储Simbiology仿真数据外，Simdata对象还提供了数据访问，绘图和分析的方法。

将变体对象的活动属性设置为true并模拟。

variantObj.Active=true；mutantData=sbiosimulate（modelObj）；

使用虚线绘制数据。另请参阅sbioplot.方便地绘制Simdata对象。

绘图（muntatdata.Time、muntatdata.Data、，“线条样式”,' - '）;传奇(mutantData。DataNames,'地点',“东北朝”)；xlabel('时间（秒）')；伊莱贝尔(“物种数量”）;网格在；

比较野生型和突变型途径中活性G蛋白物种（Ga）的行为。

gaindex = strcmp（名称，'ga'）;野生型结果的％指数[tmut，xmut]=按名称选择（修改数据，'ga'); 绘图（时间、数据（：、增益指数）、tmut、xmut、，' - ')；xlabel('时间（秒）')；伊莱贝尔(“物种数量”)；传奇({'Ga（wt）','Ga（突变）'}，'地点',“东北朝”）;网格在；

执行参数扫描

与野生型相比，突变株的g蛋白失活率要低得多(kGd = 0.004 vs . 0.11)，这解释了在上述比较中观察到的随着时间的推移g蛋白(Ga)活性较高的原因。要更详细地了解kGd的变化如何影响Ga的级别，请执行几个模拟的参数扫描，其中kGd的值在一个范围内变化。下面的示例演示了对5个kGd值的参数扫描;要增加迭代次数，请更改参数中的值邻域下面的函数。

生成五个均匀分布的kGd值，范围从0.001到0.15。

kgd值=linspace（1e-3,0.15,5）；

将参数扫描结果存储在SimData对象数组中。初始化变量以保存此数组。

Scandata = [];

为加速模拟准备模型。

sbioaccelerate（modelobj）;

循环kGd值并对每个值执行模拟。使用模型上的突变变量修改模拟期间使用的kGd值。

对于kGd=kGd值%在变体中设置所需的kGd值。variantObj.Content{1}{4}=kGd；%模拟模型，将结果存储在SimData对象中。sd=sbiosimulate（modelObj）；scanda=[scanda；sd]；结尾

斯堪的数据现在是SimData对象的五元素数组。每个对象包含参数扫描中一次运行的数据。

从SimData对象数组中提取遗传算法的时间进程，并在单个轴上绘图。下面的代码一步一步地构造情节;另外,看到sbioplot.和斯比奥苏布普洛特.

[tscan，xscan]=按名称选择（Scanda，'ga'）;fh =数字;抓住在；对于c=1:5图（tscan{c}，xscan{c}）；str=sprintf(“k=%5.3f”文本（tscan{c}（end），xscan{c}（end），str）；结尾%给情节加上注解。轴（gca（fh），“广场”)；头衔(“改变kGd值：对Ga的影响”)；xlabel('时间（秒）')；伊莱贝尔(“物种数量”）;网格在；抓住离开；

参数估计-背景

在对生物系统进行建模时，通常需要包括数值未知或仅大致已知的参数。如果系统中一个或多个物种的实验数据可用，则可以通过改变这些参数并寻找那些导致模型si之间最佳拟合的值来估计这些参数的值计算结果与实验数据吻合较好。

在本示例的这一部分中，我们将在尝试将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（表（tExpt，GaFracExpt））；data.Properties.IndependentVariableName=“tExpt”；

使用非常量参数和重复分配规则将此分数添加到模型中，而不是将此实验数据转换为Ga的绝对量。

gafracobj = modelobj.addparameter（'GAFRAC',“康斯坦特价值”, 0); GaFracRule=modelObj.addrule(‘GaFrac=Ga/（Ga+G+Gd）’,“repeatedAssignment”)

GaFracRule=SimBiology规则数组索引：规则类型：规则：1重复分配GaFrac=Ga/（Ga+G+Gd）

更改配置集上的RuntimeOptions以记录GaFrac。

configsetObj.RuntimeOptions.StatesToLog = GaFracObj;

去激活突变株。

variantobj.active = false;

模拟模型，将结果存储在Simdata对象中。

sdWild = sbiosimulate (modelObj);

获取“GAFRAC”的数据以稍后在绘图中使用。

[tWild，GaFracWild]=按名称选择（sdWild，'GAFRAC'）;

将仿真结果重新采样到实验时间向量上。

sdwildresampled =重组（sdwild，texpt，“pchip”）;

获得“Ga”物种的重新采样数据。

[~，GaFracWildResampled]=selectbyname（sdWildResampled，'GAFRAC'）;

计算r平方值，测量模拟数据和实验数据之间的拟合。

sst=norm（GaFracExpt-平均值（GaFracExpt））^2；sse=norm（GaFracExpt-GaFracWildResampled）^2；rSquare=1-sse/sst；

将模拟结果绘制于GA的实验数据。

fh=图形；绘图（tExpt，gafractexpt，“罗”)；legendText={'实验'}；头衔('适合GAFRAC'的实验数据)；xlabel('时间（秒）')；伊莱贝尔(“物种数量”）;抓住在; 情节（斜纹、野生）；legendText{end+1}=sprintf('原始，R ^ 2 =％4.2f'，rSquare）；图例（legendText{：}）；网格在；

参数估计-估计单个模型参数

从参数扫描中，我们已经看到参数kGd的值对物种Ga的时间进程有显著影响。让我们看看是否可以通过改变kGd的值来改善模型结果与实验数据的拟合。

对实验数据执行参数估计，优化KGD的值。绘图有关迭代的信息，同时优化进展，最多15个迭代。

paramToEst = estimatedInfo ('kgd'）;kgdobj = sbioselect（modelobj，“姓名”,'kgd'）;选择= OptimSet（'plotfcns'，@ OptimplotFval，“麦克斯特”15）;结果1 = sbiofit（modelobj，数据，“GaFrac = GaFracExpt”，paramToEst，...[]，'fminsearch'、选择);estvalues1 =结果1.Parameterestimates.

估计StandardError estValues1 = 1 x3表名称  _______ ________ _____________ {' kGd} 0.12142 - 0.0018542

将kGd的估计值存储在新型号变型中。

OptimVariantobj = addVariant（Modelobj，“优化kGd”)；添加内容（optimVariantObj{“参数”,'kgd',“价值”，estValues1.Estimate}）；

激活新变体并使“突变”变体失活。

optimVariantObj.Active=true；mutantVariantObj=getvariant（modelObj，'突变'）;mutantVariantObj。积极= false;

使用KGD的估计值模拟模型。

sdEst1 = sbiosimulate (modelObj);

绘制GaFrac的数据，并与之前的结果进行比较。

[t1, GaFracEst1] = selectbyname(sdEst1，'GAFRAC')；sdEst1重采样=重采样（sdEst1，tExpt，“pchip”）;[~， GaFracEst1Resampled] = selectbyname(sdEst1Resampled，'GAFRAC')；sse1=标准（GaFracExpt-GaFracEst1采样）^2；rSquare1=1-sse1/sst；图（fh）；绘图（t1，GaFracEst1，'m-')；legendText{end+1}=sprintf('kGd Changed, R^2 = %4.2f'，rsquare1）;传奇（LegendText {：}）;

从R平方值中，我们可以看出，使用新的kGd估计值，与实验数据的拟合度稍好。如果kGd的原始值只是一个粗略估计，我们可以将这些结果解释为对原始估计值的确认或对其的改进。

敏感性分析 - 背景

到目前为止，我们一直对活性G蛋白Ga的动态行为感兴趣。参数扫描显示，该物种受到控制G蛋白失活的速率常数kGd值的显著影响。通过参数估计，我们发现通过优化kGd值，我们能够更好地拟合实验结果通用航空公司的时间路线。

要问的自然问题是，模型的其他参数会影响GA水平，这些效果的大小是什么？灵敏度分析允许您通过计算相对于型号参数值或物种初始条件（“输入因子”）计算一个或多个物种（“输出”）的时间相关的衍生物来回答这些问题。

灵敏度分析-计算灵敏度

在模型中计算Ga的灵敏度。完全归一化灵敏度，以便它们可以相互比较。

停用模型上的突变变量对象，以便使用kGd在其原始值计算灵敏度。

OptimVariantobj.active = false;

在模型configset中建立灵敏度计算。

％在求解器选项中打开SensitiveSanalysis。configsetObj.SolverOptions.SensitivityAnalysis=true；%为灵敏度分析配置灵敏度输出和输入。sensitivityOpt = configsetObj.SensitivityAnalysisOptions;GaObj = sbioselect (modelObj,“类型”,“物种”,“姓名”,'ga'); 灵敏度PT.输出=GaObj；params=sbioselect（modelObj，“类型”,“参数”,“哪里”,“姓名”,“~ = ','GAFRAC'); 灵敏度点输入=参数；灵敏度归一化='满的'；

模拟模型。

sdSens=sbiosimulate（modelObj）；

从SimData对象中提取灵敏度数据，并绘制计算出的灵敏度。

[t，R，sensOutputs，sensInputs]=getsensmatrix（sdSens）；R=挤压（R）；图形；绘图（t，R）；标题(“GA关于各种参数的归一化敏感性”)；xlabel('时间（秒）')；伊莱贝尔('灵敏度'); 图例（传感器输入，'地点',“东北朝”）;网格在；

参数估计-估计多个参数

这些结果表明，Ga不仅对参数kGd敏感，而且对kGa、kRs和kRD1也敏感（其他敏感度与图上的零不可区分）。改变这四个参数可以更好地拟合实验数据。

估计这四个参数以匹配目标数据。使用先前配置的优化选项和模型中的当前参数值作为优化的起点。

选择参数kGa、kRs、kRD1和kGd进行估算。

paramsToEst=estimatedInfo({“克格勃”,“kRs”,“kRD1”,'kgd'});

如果在配置集中启用了敏感性分析选项，参数估计将忽略该选项。在求解器选项中关闭灵敏度分析以避免警告。

configsetobj.solveroptions.sensitivityAnalysis = false;结果2. sbiofit（modelobj，数据，“GaFrac = GaFracExpt”，paramsToEst，...[]，'fminsearch'、选择);estValues2 = result2。ParameterEstimates

estValues2=4x3表名估计标准错误{uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu

将四个参数的估计值存储在新模型变体中。

optimVariantObj2=addvariant（modelObj，“四个参数优化”)；添加内容（最新版本）{“参数”,“克格勃”,“价值”，estValues2.Estimate（1）}）；addcontent（optimVariantObj2{“参数”,“kRs”,“价值”，estValues2.Estimate（2）}）；addcontent（optimVariantObj2{“参数”,“kRD1”,“价值”，estValues2.Estimate（3）}）；addcontent（optimVariantObj2{“参数”,'kgd',“价值”，estvalues2.imate（4）}）;

现在用新估计的参数值模拟模型。

optimVariantObj.Active=false；optimVariantObj2.Active=true；sdEst2=sbiosimulate（modelObj）；

与以前的结果进行比较。

[t2, GaFracEst2] = sdEst2，'GAFRAC')；sdEst2Resampled=重新采样（sdEst2，tExpt，“pchip”)；[~，GaFracEst2Resampled]=按名称选择（sdEst2Resampled，'GAFRAC')；sse2=标准（GaFracExpt-GaFracEst2采样）^2；rSquare2=1-sse2/sst；图（fh）；绘图（t2，GaFracEst2，“g-”)；legendText{end+1}=sprintf(“4个常量已更改，R^2=%4.2f”，rSquare2）；图例（legendText{:}）；

由于参数估计可以自由改变四个参数，与实验数据的拟合度进一步提高。显示的优化迭代表明目标函数减小，R平方值增大。

注意，这里执行的四参数估计可能与生物学相关，也可能不相关，仅用于说明目的。

还要注意的是，将估计结果存储在变量中可以方便地在模拟不同版本的模型之间来回切换。此时有四个版本：原始版本、变体版本和基于参数估计结果的两个版本。

结论

该示例介绍了使用G蛋白信号模型进行模型建筑，模拟和分析的偶像功能的各个方面。

素质文档

万博1manbetx

试试MATLAB、Sim万博1manbetxulink和其他产品s manbetx 845

立即获得审判