文章阅读目录大纲
今天在这里和大家聊一下《生物化学系统的计算分析》这本书，生物化学系统的计算分析这本书可以说得上是我的生物信息学学习生涯的启蒙书。
大家可以上Google Books预览这本书。虽然我从后面不断地学习中了解到，其实这本书说的都是系统生物学相关的知识。但是这仍然不妨碍我为了实现相关的系统生物学模拟计算而补充学习其他相关的生物信息学知识。
系统生物学就是对生物系统的计算与数学建模过程。
我到现在任然还记得非常清楚，下面在书中的example2例子开启了我对生命网络系统的仿真模拟的好奇心；并引导了我接下来的，为了通过程序进行自动化构建下面的网络，而不断地补充其他相关的生物信息学知识（例如代谢途径注释，调控网络构建，等）的学习经历，从而诞生了GCModeller这个项目。
虽然这本书是20多年前出版的一本书，里面的知识相比于目前最新的科研进展可能显得比较陈旧了。但是这本书所介绍的S-system算法我认为是进入系统生物学领域，学习数学建模方面非常经典的一种思想方法。从S-system方程组为入口点来学习生物化学系统的仿真建模的思想，我认为还是非常值得去学习的。
为了可以让大家对生物化学系统的仿真建模有一个最基本的认识，我们现在来对于上面的截图中的途径，下面我们来具体的学习改怎么转换为S-system的动力学积分方程吧。
一般而言在S-system之中，一个变量的浓度变化的动力学过程一般由两部分组成：产生部分以及消耗的部分，也就是一个变量的动力学过程可以有一个通用的表达形式，例如：
x = a * x1 ^ a1 * x2 ^ a2 * ... - b * x1 ^ b1 * x2 ^ b2 * ...产生的部分为a部分，消耗的部分为b部分。幂指数项ai，bj等则是相应的符号xi与符号x之间的关系。可以看到，
- 如果对应的幂指数项为零，则对应的部分为1，将不会影响乘积结果，即对应的符号与x的变化无关
- 如果对应的幂指数项大于零，则表示催化效应，因为对应的幂的结果总是大于1的
- 反之如果对应的幂指数项小于零，则表示抑制效应。一个符号越大，其负的幂指数结果会越接近于零，整个乘积项将会接近于零，故表现为抑制效应。
所以对于上面的截图中的Figure 2.6的网络示意图，我们可以列出方程：
- x4为净流入，其变化量为一个常数值
- x3同样在系统内为净流入，所以其变化量也是一个常数值
- 对于x2而言，其从x4转换而来，并且转换为x1，所以a部分为x4相关，本部分为其自身降解部分。在x4转换为x2的过程，受x1催化，x3抑制，所以x1的幂指数大于零，x3的幂指数小于零
- 对于x1而言，其来自于x2的转化，并且不受其他因素影响，所以a部分只与x2相关。在消耗部分，因为受到x3的抑制，所以b部分的表达式会与x3相关，并且x3的幂指数小于零
从上面的模型可以看出，x2的消耗的表达式部分与x1的生成的表达式部分是一样的，这个是因为我们的模型要遵守自然界之中最基本的物质守恒定律。
在系统生物学里面，模拟一个动力学系统，实际上就是一个微积分方程组的求解过程。关于使用R#脚本进行动力学系统的模拟计算的原理，可以阅读之前的一篇文章：《R#语言使用龙格库塔法求解劳伦兹混沌系统》。
在这里我们仍然是使用lambda函数来表示一个微分方程组，通过这个微分方程组来构建出一个动力学系统：
[
x1 -> 2 * x2 - 1.2 * (x1 ^ 0.5) * (x3 ^ (-1)),
x2 -> 2 * (x1 ^ 0.1) * (x3 ^ (-1)) * (x4 ^ 0.5) - (2 * x2),
x3 -> 0.5,
x4 -> 1
]
|> deSolve(
y0 = list(
x1 = 2,
x2 = 0.1,
x3 = 0.5,
x4 = 1
),
a = 0,
b = 10,
resolution = 20000
)
|> as.data.frame
|> write.csv(file = `${dirname(@script)}/example2.csv`)
;执行上面的脚本求解微分方程组的积分结果，得到的曲线的值和书中的example 2示例的结果相比，比较尴尬，有点不一样。不太清楚是不是积分计算的方法差异或者是计算的分辨率的不同导致的问题：
虽然最终的结果值不太一样，但是变化趋势是与书中描述的是一样的：
X1变量首先会因为X3浓度较低，抑制力作用较低而在初期有一个较大的消耗量。在曲线上首先展现出一个被消耗的趋势，之后由于X3浓度不断累积，X1的消耗反应接近停止，X1表现出在后期不断地累积；
X2变量而言，由于其通过X4转化的过程受X3抑制，因为X3后期不断累积，所以X4转化为X2的过程在不断的受抑制，X2的转化效率低于其转化为X1的时候，X2的浓度将会不断降低
GCModeller之中内置的PLAS计算引擎
在GCModeller的simulators模块之中，定义了一个S.system
在这个程序包模块之中，提供了有一些与生物化学系统建模脚本化编程相关的api：
- snapshot 函数用来打开一个用来进行本地数据存储的文件设备，用来保存模拟器输出的每一帧符号数据，可以从下面的代码看到，这个函数接受一个符号列表用来选择性的输出结果数据
- kernel 函数用来创建一个基于S-system算法的生物化学系统动力学模拟器对象
- environment 函数就是用来配置kernel中的一些符号常量等信息
- s.system 函数就是用来向模拟器核心加载我们所定义的S-system系统方程
- run 最后run函数就是执行整个模拟器的运行
PLAS模拟器脚本示意
下面我展示了demo代码，查看如何通过GCModeller所提供的PLAS系统相关的api来进行动力学系统模拟计算的脚本化编程：
imports "S.system" from "simulators";
# set y0 of the ODEs system
const symbols = list(
x1 = -100,
x2 = -100000,
x3 = 0,
x4 = 10
);
using data.driver as snapshot("./atkinson.csv", symbols = names(symbols)) {
data.driver
|> kernel(S.script("Atkinson system"))
|> environment(
beta1 = 30,
beta3 = 30,
beta4 = 1,
lamda1 = 2,
lamda3 = 2,
alpha1 = 20,
alpha2 = 20,
alpha3 = 1,
a = 1,
n1 = 4,
n2 = 5,
n3 = 1,
is.const = TRUE
)
|> environment(symbols)
|> s.system([
x1 -> beta1 * (lamda1 * (1 + alpha1 * (x4 ^ n1) / (1 + x4 ^ n1)) - x1),
x2 -> x1 - x2,
x3 -> beta3 * (lamda3 * (1 + alpha2 * ((x4 / a) ^ n2) / (1 + (x4 / a) ^ n2)) * (1 / (1 + x2 ^ n3)) - x3),
x4 -> beta4 * (x3 - x4)
])
|> run(ticks = 4.5, resolution = 0.1)
;
}可以看到，执行上面的脚本代码，上面的S-system动力学方程组输出了这样的一个无限震荡并且信号放大的动力学系统。
- 【MZKit】简单自动化组织分区 - 2023年11月5日
- 【MZKit教程】质谱成像原始数据文件查看 - 2023年6月29日
- 生物序列图嵌入算法 - 2023年6月29日
No responses yet