一种用可逆单分子反应实现马尔科夫链的设计方法与流程

本发明属于化学反应网络(crns)计算领域，尤其涉及一种用可逆单分子反应实现马尔科夫链的设计方法。

背景技术：

通常情况下，马尔科夫链被用来分析化学反应网络中的随机过程。但是，s.a.salehi,m.d.riedel,andk.k.parhi在发表于ieeeinternationalconferenceondigitalsignalprocessing(dsp),2015,pp.689–693的论文“markovchaincomputationsusingmolecularreactions”，以及m.cardona,m.colomer,j.conde,j.miret,j.miro,anda.zaragoza，在发表于biosystems,vol.81,no.3,pp.261–266,2005的论文“markovchains:computinglimitexistenceandapproximationswithdna”中做了相反的工作：用dna反应来估算马尔科夫链的稳态分布。s.a.salehi等人首先设计了形式化学反应网络，m.cardona等人将dna链被用来表示马尔科夫链的端点和边，但是上述的方法都没有实现连续时间马尔科夫链的计算。所以本发明致力于设计化学反应网络(crns)来计算包括连续马尔科夫链的计算，并在一定程度上降低了前人的复杂度。

基于常微分方程组(odes)的确定性模型能够很好地表达化学反应网络的动力学特性，所以本发明使用ode模型来进行仿真。根据质量作用定律，反应发生的速率正比于反应物的浓度和速率常数。比如对于反应a+b→c+d，ode模型如下：

2010年，d.soloveichik等人提出，任何双分子或单分子的反应可以映射到dna链置换反应，所以可以随意设计假想的化学反应网络。值得注意的一点是量化的系统下的速率常数和浓度的度量保持不变。所以，本发明中的仿真浓度和时间都是没有单位的。最后一点，本发明中所有反应的速率常数都是1，在下文中不再赘述。

技术实现要素：

发明目的：针对现有技术中用dna反应估算马尔科夫链稳态分布的方法不能实现对连续时间马尔科夫链的计算，本发明提出了一种用单分子可逆反应组成的化学反应网络实现马尔科夫链的设计方法，该方法同时适用于离散时间的马尔科夫链和连续时间的马尔科夫链。

技术方案：为了实现上述目的，本发明中用单分子可逆反应实现马尔科夫链的设计方法，包括以下步骤：

(1)设计化学反应网络，用不同的反应物种类代表目标马尔科夫链中的不同状态，根据目标马尔科夫链各状态的初始概率分布设置相应反应物的初始浓度；

(2)用反应物之间的反应速率常数表示转移概率，根据转移概率设置反应速率常数的值；

(3)用可逆单分子反应表示状态转换；

(4)用设计好的反应网络计算马尔科夫链的稳态分布，所有反应物的最终浓度即为稳态分布。

其中，步骤(1)中反应物的初始浓度设置必须等于或正比于目标马尔科夫链中相应状态的初始概率分布。

其中，步骤(2)中所有反应物之间的反应速率常数值必须等于或正比于目标马尔科夫链中相应状态之间的转移概率。

有益效果：本发明中用单分子可逆反应实现马尔科夫链的设计方法通过构建化单分子可逆反应组成的化学反应网络与马尔科夫链之间的关系，利用单分子反应物的反应速率常数表征两两状态之间的转移概率，有效地利用了化学反应可逆的本质，有效地减少了所需反应物的种类和反应个数。

附图说明

图1为本发明中用单分子可逆反应实现马尔科夫链的设计方法的流程图；

图2为一个离散马尔科夫链及实现该链的化学反应网络，图2(a)为一个离散马尔科夫链的状态转换图，图2(b)为实现图2(a)中马尔科夫链的化学反应网络；

图3为图2中马尔科夫链的ode仿真结果；

图4为一个连续马尔科夫链及实现该链的化学反应网络，图4(a)为一个连续马尔科夫链，图4(b)为实现图4(a)中马尔科夫链的化学反应网络；

图5为λ＝0.5时图4中马尔科夫链的仿真结果；

图6为λ＝1时图4中马尔科夫链的仿真结果。

具体实施方式

如图1所示，本发明中用单分子可逆反应实现马尔科夫链的设计方法，包括以下步骤：

(1)用不同的反应物种类νi(i＝1,2,...,k)νi,i＝1,2,…k代表不同的状态，根据目标马尔科夫链(待求解稳态分布的马尔科夫链)的初始概率分布设置反应物的初始浓度；即所有反应物的初始浓度必须等于或正比于目标马尔科夫链的初始状态分布；

(2)用反应速率常数表示转移概率kij(i＝1,2,…,kj＝1,2,…,k)，根据转移概率设置反应速率常数的值；

(3)用可逆单分子反应表示状态转换；

(4)用设计好的反应计算马尔科夫链的稳态分布，所有反应物的最终浓度即为稳态分布。

上述步骤(4)中，对于连续马尔科夫链，反应物浓度随时间的变化曲线即为该链的瞬时解，所有反应物的最终浓度即为稳态分布。

本发明方法以反应物最终的浓度表示了稳态概率分布，在连续马尔科夫链中，任意时刻t反应物的浓度表示了t时刻所有状态的概率分布。也就是说，该方法能够用化学反应网络计算马尔科夫链的稳态分布和瞬时解。

针对离散时间马尔科夫链(dtmc)，本发明方法具体介绍如下：

离散时间马尔科夫链(dtmc)是对所有n∈ν0,si∈s都满足

p(xn+1＝sn+1|xn＝sn,xn-1＝sn-1,...,x0＝s0)＝p(xn+1＝sn+1|xn＝sn)的随机过程{x0,x1,...,xn+1,...}。在齐次的情况下，从状态i到状态j的转移概率定义为pij＝p(xn+1＝j|xn＝i)。向量ν(n)＝(ν0(n),ν1(n),ν2(n),...)表示n时刻所有状态的概率分布。

为了实现相应的功能，该设计方法是用不同的化学反应物种类模拟每一种状态。每一个单分子反应表示一种状态转换。所有的速率常数值都等于或正比于转换概率。所有的反应物初始浓度都等于或正比于目标马尔科夫链的初始概率分布。dtmc和crn的映射关系总结在表1。实现计算的反应网络如下：

表格1dtmc和crns的简单映射关系

针对连续时间马尔科夫链(dtmc)，本发明方法具体介绍如下：

连续时间马尔科夫链(ctmc)是对任意

ti∈r0⁺,0＝t0＜t1＜…＜tn＜tn+1,si∈s满足

的随机过程{xt:t∈t}。向量π(u)＝(π0(u),π1(u),π2(u),...)表示u时刻所有状态的概率分布。实时转移速率定义为

为了实现目标功能，连续马尔科夫链的设计方法与离散情况时类似。ctmc和crn的映射关系总结在表2。不同的化学反应物种类模拟每一种状态。每一个单分子反应表示一种状态转换。所有的速率常数值都等于或正比于实时转移速率。所有的反应物初始浓度都等于或正比于目标马尔科夫链的初始概率分布。实现计算的反应网络如下：

表格2ctmc和crns的简单映射关系

图2(a)是马尔科夫链的状态转换图，从图中可以看出，这条链一共有100个状态，分别是状态1,2,3,…100。如果当前是状态1，那么下一个时刻状态不变。如果当前是状态2，那么下一个时刻有0.6的概率跳转到状态1，有0.4的可能性跳转到状态3，以此类推。现在利用上述的映射方法用化学反应网络计算该链的稳态分布，用100种反应物ν1,ν2,…，ν100表示该链的100种状态。假定初始状态是99，所以ν99的初始浓度是1，其余物质的初始浓度均为0。图2(b)为实现计算的反应网络，根据转移概率设置反应速率常数的值，本实施例中设置反应速率常数等于对应的转移概率。

本发明中用了确定性的化学反应网络模型来仿真，以证实设计方法的准确性。具体来说就是用常微分方程组(odes)进行仿真。反应网络设计完毕后，根据物质作用定律，列出相应的ode方程，解出所有物质浓度随时间的变化。在反应网络设计完毕后，使用来源于http://users.ece.utexas.edu/～soloveichik/的仿真软件，crnsimulator，通过mathematica进行仿真。输入所有的反应和反应物初始浓度，就可以得到仿真图像，如图3。从仿真的结果来看，ν1和ν100的最终浓度分别趋于0.33和0.67，其余物质的最终浓度为0，故在图像中省略。说明这条链的稳态分布包括ν1和ν100，概率分别为0.33和0.67，其余状态为0。这与数学上解出的答案相符。

图4(a)中的马尔科夫链有无数个状态，其中每个状态都以相同的速率λ跳转到下一个状态。由于实际中可以利用的化学反应物有限，所以用六个反应物πi,i＝0,1,2,…5近似表示这条链，从连续马尔科夫链及其瞬时解的定义来看，此时πi,i＝0,1,2,3,4五个状态的瞬时解与有无数个状态时的前5个状态的解完全一样，所以仿真时只显示这五个状态的图像。假定初始状态为π0，所以π0的初始浓度是1，其余物质的初始浓度是0。图4(b)中设计的反应网络，其反应速率与转移速率相同，均为λ。仿真过程与上述离散的情况一样。图5展示了当λ＝0.5时图4中马尔科夫链的仿真结果，图6展示了当λ＝1时图4中马尔科夫链的仿真结果，其中t时刻所有物质的浓度即为该链的瞬时状态分布。从图像上观察，仿真结果与通过数学方法解除的瞬时解一样。

以上仅是本发明的优选实施方式，应当指出以上实施列对本发明不构成限定，相关工作人员在不偏离本发明技术思想的范围内，所进行的多样变化和修改，均落在本发明的保护范围内。