柱坐标系下偏微分方程与浓度模型的模拟及应用

理论计算机科学电子笔记268（2010）49-59www.elsevier.com/locate/entcs柱坐标系下偏微分方程和CTMCs与浓度Andrea Degashte1和Mu Zarey Calder2格拉斯哥大学计算机科学系英国格拉斯哥摘要本文提出了一种在圆柱坐标系下，由Fick定律的偏微分方程。由此产生的模型抽象的摩尔浓度，离散的水平，和空间位置，离散区室。 我们将结果应用于一氧化氮在人体血管中的作用，并用PRISM工具中的模拟进行说明关键词：几何空间，反应微分方程，柱坐标，偏微分方程，常微分方程，连续时间马尔可夫链1引言诸如过程代数和其他演算[15，19，11，6]、重写规则[12，5]或语言[8，18]等形式主义可以用于改进生化反应系统的建模和分析。通常，用这些方法定义的模型使用数学技术，如普通微分方程（ODE），连续时间马尔可夫链（CTMC），CTMC与水平[9]或蒙特卡洛模拟作为底层的具体语义。在某些情况下，可以从相同的形式主义中导出多个数学语义，例如Bio-PEPA [11]，并且它们可以彼此相关，以便对结果进行更鲁棒的解释[9]。最近，越来越多的兴趣已经得到了这样的形式主义在空间中的位置和运动的整合。空间位置和生物化学物种的多样性可以用多种方式来表示。例如，空间可以是拓扑的，即分层位置，或者是几何的，即空间的坐标系。1电子邮件地址：andrea@dcs.gla.ac.uk2电子邮件：muffy@dcs.gla.ac.uk1571-0661 © 2010 Elsevier B. V.在CC BY-NC-ND许可下开放访问。doi：10.1016/j.entcs.2010.12.00550A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）49空间位置[13]。分子的离散性可以在微观水平上用随机游走来描述，或者在宏观水平上用菲克离散性定律来描述在这种情况下，数学模型通常是高度具体化的，简化了所用方程组的假设。这方面的一个例子是一氧化氮（NO）运输和血管中可用性的各种数学模型（完整综述见[20]）。这些模型都使用偏微分方程组（PDE），在一维柱坐标系中用Fick定律表示微分。更详细地，血管被建模为具有表示组织层的同心隔室的圆柱体。分析的主题通常是NO可以扩散的距离，从产生NO的组织层。为了计算这样的距离，可以简化模型，只考虑一个维度，半径，因为其他维度是不变的。我们的目标是得到一个近似模型的NO运输和可用性在血管中的CTMC水平。具有水平的CTMC是其状态由以离散水平表示的每种物质的浓度表征的CTMC。这样做的动机是可用的额外分析，例如，在不同随机程度下的稳健性测试。此外，在这种语义方面的偏微分方程的近似值确保了我们可以使用Bio-PEPA，它是基于CTMC的水平，以及为它开发的工具来建模和分析这种情况。在我们这样做之前，拼图中缺少的部分是一维柱坐标中的离散度在CTMC中的推导。在本文中，我们提出了这样的推导，其主要新颖之处是：• 所得到的CTMC与水平的比率直接从使用生物化学物质的二聚体的二聚体常数导出 这种推导，在情况下微不足道，在圆柱坐标的情况下，需要额外的假设;• CTMC与水平的比率不仅取决于物种的浓度和隔室的体积，而且还取决于空间位置。本文的组织结构如下。在第2节中，我们提出了推导的偏微分方程的CTMC水平。在第3节中，我们给出了一个例子。在第4节中，我们提到了相关的工作，而结论和未来的工作在第5节。2一维圆柱形容器中的扩散：偏微分方程、常微分方程和CTMC模型本节中的推导是指圆柱坐标中的一维模型。虽然考虑的唯一尺寸是半径，但必须注意，距离中心r的每个点的浓度代表半径为r的圆的圆周上每个点的浓度。圆柱的长度可以忽略，因为它是我们推导的不变量。因此，我们把面积和体积当作同义词使用在我们给出推导的细节之前，我们给出一个简短的概述。首先，我们介绍了偏微分方程，我们推导出一个数值逼近的常微分方程。这A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）495122通过将空间分割成段并通过相对于径向位置的一阶和二阶导数的数值近似来获得其次，我们描述的微分方程与隔间，其中的传输反应的速度 第三，我们证明，如果我们选择适当的动力学常数的传输质量作用定律，所得到的方程是相同的偏微分方程的近似。第四，我们从常微分方程中推导出具有水平的CTMC，类似于[10]中已经提出的推导。2.1偏微分方程物种S的建模分歧在宏观层面上由菲克方程定义，这是一个偏微分方程（PDE），其形式为δ[S] =DSδt中国（1）其中DS是物质S的扩散系数，[S]是S（单位：摩尔，M）和ε2是拉普拉斯算子，根据坐标系的不同，可以以不同的方式解释。由于我们考虑了局部生物化学相互作用的扩散，我们使用以下反应扩散方程[16，4]：δ[S] =DSδt中文（简体）其中React表示涉及物质S的其他反应。我们假设圆柱坐标，浓度仅随径向位置变化，并且仅由于扩散或生化反应。因此，方程（2）可以简化为圆柱坐标中的一维形式，其中唯一考虑的维度是半径：δ[S]（r）. δ2 [S] 1δ[S]δt=DS·边界条件是：δr2 +r·δr（3）δ[S]=零δrr=0δ[S]δrr=R= 0（4）其中r= 0表示坐标系的中心，而R是所考虑的圆形区域的半径。在时间上的每个时刻t，我们可以从初始浓度分布f（r）开始计算S在半径上任何点r处的浓度。现在考虑如何数值求解方程（3）。首先，我们将半径划分为长度为Δr =R/K的K段。 每个段i都与变量[Si]相关，i = 1，.，K，代表该区域的平均浓度。 其次，我们使用[Si]计算[S]在径向位置的一阶和二阶导数的近似值。这些近似表示第i段中点的导数，距离坐标系中心r= Δr（2i− 1）/2。使用中心有限差分计算导数52A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）49Σ·±δt（Δr）2Δr（2i− 1）/2（Δr）2（2i−1）一期+1我（2i−1）i−1δt（Δr）2（2i−1）我（2i−1）KK=K−1K=K K−1δt（Δr）221（八）方法：δ[S]（r）δrδ[Si]=δr[Si+1]−[Si−1]2Δrδ2[S]（r）δr2δ2 [Si]δr2为[Si+1]− 2[Si]+ [Si−1]（5）（Δr）2我们现在可以使用等式（5）中的近似来重写等式（3δ[S]（r）δ[Si]=D·[Si+1]−2[Si] +[Si−1]+1[Si+1]−[Si−1]React2Δr=DS·。. 1+1[S]−2[S ]+。1−1[S]±React重写最后一个方程，我们得到最终的数值近似：δ[Si]=DS·。2i[S]−2[S ]+（2i−2）[S]Reacti= 2，.， （K − 1）为了也针对i= 1和i=K写出等式（6），我们需要采用边界条件（等式（4）），从该边界条件我们获得：[S1]−[S0]= 0[SK+1]−[SK]= 0ΔrΔr因此，等式（5）中的近似变为：δ[S1]=δr[S2]−[S1]2Δrδ2[S1]δr2= [S2]−[S1]（Δr）2δ[S] ]−[S]δ2[S][S]−[S]（七）δr2Δrδr2（Δr）2利用等式（7），我们推导出两个附加等式：δ[S1]=2DS·。[S]−[S]±Reactδ[SK]=2K−2·DS·。[S]−[S]React作为最后一步，我们使用初始条件f（r），r ∈ [0，R]导出时间t = 0时[Si]的值。 注意半径为r的圆周上的每一点都有浓度f（r）。这意味着当我们计算[Si]时，第i段中的每个点都有不同的权重，即该段的平均浓度由于Vi = π（Δ r）2（2i−1）是环的面积，其平均浓度由[Si]表示，因此t = 0时的平均浓度由下式给出：1[Si]（t= 0）=我Δr（i）Δr（i−1）∫VδtS一期+1i−1（六）δt2K− 1（Δr）2K−1KA. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）49532πr·f（r）dr（9）54A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）491δt21我 δt一期+1我我i−1KδtK−1KKk（i-2），（i-1）k（i-1），iRki，（i+1）k（i+1），（i+2）S（i-2）S（i-1）SiS（i+1）S（i+2）k（i-1），（i-2）ki，（i-1）k（i+1），ik（i+2），（i+1）i-2 i-1 i i +1 i +2图1.一、在隔间中划分空间2.2常微分方程在本节中，我们根据常微分方程（ODE）推导出方程（3）的近似，其中我们将空间离散为具有体积Vi（i = 1，...，K）。 为了表示存在隔室的物种S，我们需要为每个隔室Ci指定一个变量Si。我们用[Si]来确定Ci中物质S的平均浓度浓度[S i]可以从一个隔室C i迁移到相邻的隔室，即C i−1或C i+1，并且迁移以质量作用定律给出的速度发生，动力学常数为k i，i−1（i = 2，...， K）或k i，i+1（i = 1，.， （K− 1））（见图1）。特别地，我们证明了动力学常数可以从微分方程常数DS导出，并且对于给定的K，偏微分方程和常微分方程的数值解是等价的。上述ODE系统由以下等式组成：δ[S1]V·=k[S ]−k[S ]±V·Reactδ[Si]V·=k[S]−k[S]−k[S ]+k[S]±V· 反应i= 2，.， （K − 1）δ[SK]V·=k[S]−k[S ]±V·反应其中体积Vi=π（Δr）2（2i− 1）。然后我们可以重新排列上述方程：δ[S1]=k2， 1[S ]−k1， 2[S ]± React（10）δt V12V11δ[Si]=ki+1，i[S]−ki，i+1[S]−ki，i−1[S ]+ki−1，i[S]±Reactδt Vi一期+1ViiViiVi ii−1（十一）δ[SK]=kK−1，K[S]−kK，K−1[Si= 2，.， （K − 1）]± React（12）δtVKK−1VKK在这一点上，我们选择动力学常数，参数在方程（10），（11）和（12）中代入，得到方程（6）和（8）。我们推导出这些二，一一、二1i+1，ii，i+1i，i−1i−1，i我K−1，KK，K−1A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）4955（Δr）2一期+1（2i−1）我（2i−1）（2i−1）δt（Δr）2（2i−1）我（2i−1）δt（Δr）221δt2K− 1（Δr）2根据等式（6）的检查得到的常数：ki，i+1S2i= 2iπDi（Δr）2（2i− 1）Si= 1，.， （K − 1）（13）VDS （2i−2）=（2i− 2）πDi= 2，.，中文（简体）i，i−1i（Δr）2（2i− 1）S注意，ki，i+1=ki+1，i。然后，我们可以将等式（13）和（14）替换为等式（11）：δ[Si]=（2（i+1）−2）πDS[S]−2iπDS[S]我δt π（Δr）2（2i−1）一期+1π（Δr）2（2i− 1）- -一种 （2 i− 2）πDS[S ]+ （2（i−1））πDS[S我]±Reactπ（Δr）2（2i−1）π（Δr）2（2i−1）=DS·2 i[S] −（2 i +2 i − 2）[S]+（2 i −2）[S]React最后再重新安排一下：δ[Si]=DS·。2i[S]− 2[S ]+（2i− 2）[S]Reacti= 2，.， （K − 1）这与等式（6）相同。以类似的方式，我们可以从等式（10）和（12）开始导出两个附加等式δ[S1]=2DS·。[S]−[S]±Reactδ[SK]=2K−2·DS·。[S]−[S]React其与等式（8）相同。初始条件的推导完全如偏微分方程所示因此，我们已经表明，我们可以推导出常微分方程的隔间从偏微分方程在柱坐标。这是可能的，因为我们已经发现在偏微分方程中的扩散常数DS和在常微分方程中的质量作用传递反应的动力学常数ki，j之间的对应关系在这一点上，我们已经从[10]中知道，我们可以导出相应的CTMC与水平。2.3具有集中在前两节中，我们将连续空间偏微分方程模型与离散空间常微分方程模型联系现在我们考虑进一步的离散化：我们将后者的连续浓度与水平模型的CTMC的离散浓度i−1一期+1i−1K−1K56A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）49联系起来[9，10]。该模型的组织与刚才介绍的模型相似：空间被划分为隔间Ci，体积为Vi=π（Δr）2（2i− 1），每个隔间代表一个环，其中物质S的平均浓度由[Si]A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）4957我 δt我δt我一期+1δtΔt我（i = 1，...， K）。然 而 ，该浓度并不像在ODE模型中那样以连续形式表示，而是通过离散水平Δ S iΔ = 0，.， N i，其中N i为最大水平。浓度的量可以在任何时候使用关系式[Si]= Si·hi进行评估，其中hi称为步长，是由一个水平表示的浓度的量。现在我们将定义CTMC的状态和跃迁，其能级来自于前几节中给出的模型。具有水平的CTMC的状态由水平向量σ =（Σ S 1，.，S K为了使CTMC的状态空间是有限的，对于每个变量Si，最大浓度Mi是固定的，以Ni个间隔划分，表示hi摩尔浓度，其中hi= Mi/N i，i = 1，.， K.利用生化反应及其速度，从常微分方程模型中推导出CTMC随能级的跃迁及其随时间的活化。我们使用以下附加符号：R i，i+1：S i→ S i+1v i，i+1= k i，i+1 [S i]i = 1，.，（K − 1）Ri，i−1：Si → S i−1vi，i−1 = ki，i−1 [S i]i = 2，.，K其中反应Ri，j表示Si转化为Sj（即S从Ci迁移到Cj），而vi，j是反应Ri，j的速度，以摩尔/s表示，j∈ {i+ 1，i− 1}。此外，我们假设当反应R i，i+1或R i，i-1发生时，CTMC将从状态σ=（？S1？，.，Si−1状态σJ=（S1，.，SiSKSi−1<$+1，SK如果<$Si<$= 0或<$Sj<$= N j，则反应R i，j不能发生。现在考虑单个反应Ri，i+1的常微分方程。它由两个互补方程组成：V·δ[Si]=−k[S]Vδ[Si+1]·=k[S]i = 1，.， （K − 1）（15）此外，等式（15）可以写成以下差分形式：δ[Si]V·V· Δ Si·hi=−k·S·hiδtiΔti，i+1ii（十六）δ[Si+1]V·V·Δ Si+1·hi+1 =k·S·h我们对CTMC态跃迁的假设是：当反应Ri，i+1发生时，Si的一个能级被消耗，Si+1的一个能级被产生.这意味着在等式（16）中，Δ<$ Si<$ =− 1且Δ<$ Si+1<$= 1。 现在请注意，方程（16）中唯一的未知项是Δ t，它可以被视为将电平S i转换为电平Si+1所需的平均时间。 所以我们有：Δt =Vi·hi=Vi+1·hi+1i= 1，.， （K − 1）（17）ki，i+1·Si·hiki，i+1· Si·hii，i+1一期+1i，i+1一期+1i，i+1我58A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）49重新排列等式（17），我们得到以下等式：h=hi+1·Vi+1=h（2i +1）·i= 1，.， （K − 1）IVi一期+1（2i− 1）作为一个主要的结果，我们有，一旦一个步长hi被选中，其他的步长是自动派生的我们建议先计算h1，虽然也可以从其他的隔间开始特别地，从C1开始将确保所有其他的h i（i = 2，.，K）将小于h1。然后我们可以使用ri，i+1= 1/Δ t作为反应Ri，i+1所需时间的指数分布的参数。 通过对方程（17）的处理，我们得到：ri，i+1 =ki，i+1SiViki，i−1<$Si<$i= 1，.， （K −1）（十八）ri，i−1=Vi= 2，.，K2.4附注事项虽然我们已经证明了偏微分方程的数值解和微分方程的常微分方程模型是等价的，但还需要进一步考虑在PDE的情况下，K将对建模者隐藏，并且通常相当大，以便获得尽可能接近分析解的输出。当转换为常微分方程时，建议使用较低的K，因为建模者必须直接处理隔间，并且其中大量的隔间很难管理。从常微分方程到具有能级的CTMC，我们注意到CTMC的状态空间取决于K，另外还取决于最大能级数Ni。直观地，较大的K和Ni产生马尔可夫链，其输出越来越趋向于原始PDE的输出最后，我们注意到CTMC的复杂性，再现PDE输出的能力，以及由于浓度离散化的随机效应之间的紧张关系。管理这种紧张关系是建模者的工作3例如我们现在将注意力转向一个例子，其灵感来自一系列关于人体血管中一氧化氮（NO）代谢的出版物特别地，在[17]中定义了半径R为138μm的血管;我们认为NO的扩散常数DNO= 3300μm2s−1。作为t = 0时的初始浓度函数f（x），我们选择：Γ（α+β）f（r，α，β）=Γ（α）·Γ（β）·. rα−1R.R−r<$β−1R·10其中α= 1，β= 3，定义为[0，R]，单位为μ摩尔。这种选择是基于我们在文献中的NO浓度的经验，作为测量的结果或在其他数学模型中观察到它在α中是参数的我·A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）4959ODE分组演化3025201510500 0.5 1 1.5 2 2.5 3时间s图二.微分器前3秒的PDE输出（左）和ODE输出（右）。线是沿半径的样品点处的浓度（左）和隔室Ci（i= 1，...，10）（右）。和β，以允许从相同的函数开始生成全范围的初始条件。注意，到目前为止的信息足以在一维柱坐标中求解方程（1我们使用模拟器FlexPDE [2]（图2左侧）。通过将隔室的数量定义为K= 10，我们也可以计算ODE解。对于这项任务，我们使用了模拟器Copasi [1]（图2右侧，其中绘制了每个隔室的平均浓度CTMC模型的实施也需要很少的额外信息。我们定义最大能级数N1 =10，最大浓度Mi =40 μmolar，i = 1，.，K.这里我们使用了Prism模型检查器[3]。结果是CTMC具有1. 13个国家和2。4· 1014个转换。 大型CTMC实际上是预期的，因为状态的数量相对于模型中物种的数量呈指数增长。为了分析瞬态特性，链，随机模拟或状态空间减少技术的建议。作为对链属性的首次探索，我们使用了随机模拟，从100次运行中获取模型输出的平均值和标准差。一些模拟结果如图3所示。4相关工作从ODE到CTMC的转换与浓度水平的关系[7]的文件。这在[9]中得到了进一步的研究，其中引入了两种方法之间更坚实的理论联系。在[10]中最终考虑了隔室之间的运输。我们将微分方程从偏微分方程转化为带隔室的常微分方程的出发点是[14]，它考虑了反应微分过程的随机模拟然而，在[14]中只考虑笛卡尔坐标。浓度梯度60A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）49[S1][S5][S10][S5] ODE浓度梯度浓度梯度单次模拟，S1，10个级别4035100次模拟，S1有10个级别403530 3025 2520 2015 1510 105 500 0.5 1 1.5 2 2.53时间s00 0.5 1 1.5 2 2.5 3时间s100个模拟，S5有90个级别4035100个模拟，S10有190个级别403530 3025 2520 2015 1510 105 500 0.5 1 1.5 2 2.53时间s00 0.5 1 1.5 2 2.5 3时间s图三. CTMC随机模拟，100次运行的平均值和标准差。S1、S5和S10是隔室C1、C5和C10的平均浓度。显示了100次运行的标准差5结论和今后的工作本文从一维圆柱坐标系中的偏微分方程出发，导出了含能级CTMC的Fick微分定律。作为中间步骤，我们将偏微分方程转换为带隔室的普通微分方程，其中输运速度用质量作用动力学定律实现。这种推导的新颖之处在于动力学常数直接从摩擦学常数导出，并且它们依赖于径向位置。虽然这种推导在笛卡尔坐标的情况下是微不足道的，但在圆柱坐标的情况下必须考虑额外的假设。最后通过一个例子说明了结果的一致性，并给出了偏微分方程、常微分方程和CTMC的在未来，我们计划开发和分析一个完整的CTMC与水平模型的NO运输和生物利用度在血管中，其中的差异是使用这里提出的推导实现确认Andrea Degashte获得格拉斯哥大学的开尔文勋爵/亚当·史密斯奖学金和EPSRC资助的SIGNAL项目的支持。[S1] ODE[S10] ODE浓度梯度浓度梯度A. 德加斯托湾Calder/Electronic Notes in Theoretical Computer Science 268（2010）4961引用[1] Copasi网站（2010年）。URLhttp://www.copasi.org/tiki-index.php[2] FlexPDE网站（2010年）。URLhttp://www.pdesolutions.com/[3] 棱镜网站（2010年）。网址http://www.prismmodelchecker.org/[4] Berg，H. C.的方法，“Random Walks in Biology,” Princeton University Press,[5] Blinov，M. L.，J. R.费德湾Goldstein和W. S. Hlavacek，BioNetGen：基于分子结构域相互作用的基于规则的信号转导建模软件，Bioinformatics20（2004），pp. 3289-3291。[6] 博尔托卢西湖Stochastic concurrent constraint programming，Electronic Notes in Theoretical ComputerScience164（2006），pp. 65比80[7] 考 尔 德 ， M. ， S. Gilmore 和 J. Hillston ， 使 用 随 机 过 程 代 数 PEPA 的 影 响 ， 计 算 机 科 学 讲 义 4230（2006），pp.1-23[8] 卡尔佐内湖，F. Fages和S. Soliman，Biocham：An environment for modeling biological systems andformalizing experimental knowledge，Bioinformatics22（2006），pp. 1805-1807年。[9] Ciocchetta，F.，A. Degasstrom，J. Hillston和M. Calder，关于CTMC和源自bio-pepa的ODE模型的一些研究，ENTCS229（2009），第100页。145-163.[10] Ciocchetta ， F. 和 M. L. Guerriero ， 模 拟 bio-pepa 中 的 生 物 隔 室 ， 电 子 笔 记 理 论 。 Comput. Sci. 227（2009），pp.七十七比九十五[11] Ciocchetta ， F. 和 J. Hillston ， Bio-PEPA ： 生 物 系 统 建 模 和 分 析 的 框 架 ， 理 论 计 算 机 科 学 410（2009），pp.3065-3084[12] 达诺斯，五，J. Feret，W.丰塔纳河Harmer和J. Krivine，细胞信号的基于规则的建模，计算机科学讲义（2007），pp.17比41[13] Degasstrom，A.和M. Calder，关于生物系统梯度扩散模型的形式化，PASTA研讨会2009（2009）。URLhttp://www.dcs.gla.ac.uk/www.example.com~andrea/files/DegasperiCalder_Bio-PASTA2009.pdf[14] 埃尔班河，J. Chapman和P. Maini，反应扩散过程的随机模拟实用指南，ArXiv e-prints（2007），pp.1-35[15] Hillston，J.，[16] Jones，D. S.和B. D.李明，[17] Lamkin-Kennard，K.一、D. G. Buerk和D. Jaron，微循环中NO和O2之间的相互作用：数学分析，微血管研究68（2004），pp.38比50[18] Pedersen，M. 和G. Plotkin，生物化学系统语言，生物信息学第5307（2008）号来文，第5307（2008）页。63比82[19] 普里阿米角和p.Quaglia，用于生物相互作用的β-结合剂，LNCS3082（2005），pp.20-33[20] Tsoukias，N.M.，微循环中的一氧化氮生物利用度：数学模型的见解，微循环15（2008），pp。813-834