Mono Visual Basic 10 .NET4.5 架构师中级开发 Windows .NET Visual Basic Linux

使用 VB 对生化系统进行建模

谢桂纲

5.00/5 (9投票s)

2013年10月6日

CPOL

9分钟阅读

29820

725

S-系统方程的数学方法，用于模拟生物化学网络系统。

引言

系统生物学是一门研究生物系统以及系统中组件之间相互作用所产生的行为的学科。系统生物学研究的主要工作是建模生物系统，然后使用程序模拟该模型。由于我正在努力进行细菌萎蔫黄单胞菌致病变种8004的调控信号转导网络分析，这是我的微生物遗传学研究工作的一部分，我需要一个程序来进行网络模拟和数据分析。因此，我重写了由Antonio E.N.Fereira 基于E.O.Voit 的研究工作编写的PLAS程序，以建立我自己的模拟系统，并在此发布我的编码工作，希望能帮助其他生物研究人员。

虽然PLAS目前没有LINUX版本，但PLAS软件使用wine程序在LINUX平台上运行良好。

PLAS软件网站：http://enzymology.fc.ul.pt/software/plas/

背景

生化反应系统和S-系统功能简介

考虑一个生化系统有N（N>=1）种参与反应网络的底物，且系统温度和体积不变。因此我们可以用向量X(t)来表示系统在任何时间t的状态

X(t)=(x₁(t), x₂(t), ..., X_N(t))

x_i(t) 表示一种底物的数量或浓度（i=1,2,…,N）。

根据E.O.Voit的工作观点，一个反应过程可以用S-系统函数表示

x_i^'= α_i∑(x_j)^gij-β_i∑(x_j)^hij

x_j 代表一种底物的浓度，参数α_i , β_i , g_ij和h_ij 代表与反应x_i'相关的参数，它们都可以从实验室实验工作中获得，而x_i' 也代表底物x_i 的生成率。

因此，一个像这样的反应系统

x₁ -> x₂ -> x₃

这样一个链式反应可以用S-系统函数表示

因为每个方程的左侧x_i' 表示底物生成率，所以在计算方程右侧的值后，将其乘以时间单位，即可得到该底物的数量或浓度变化值。然后，我们得到底物的新浓度，如下表达式所示：

x_i=x_i+x_i'*△t

因此，我们可以使用新浓度在下一个迭代循环中计算其他底物。

所有这些知识和更详细的信息都可以在E.O.Voit早年撰写的《生化系统计算分析》一书或其科学研究文章中找到。

使用代码

使用VB构建生化网络系统模拟器

生化系统中的元素

生化系统中的元素可分为三类：底物、反应和系统扰动。

底物

您可以将底物视为生化系统网络中的一个节点。它有两个主要属性来代表这个节点：名称和浓度或数量值。我们在模型中定义的底物类仅用于存储计算结果。

它的数据结构定义正如你所见，非常简单

Public Class Var
 
    <Xml.Serialization.XmlAttribute> Public Name As String
    <Xml.Serialization.XmlAttribute> Public Value As Double
    <Xml.Serialization.XmlAttribute> Public Title As String
    <Xml.Serialization.XmlAttribute> Public Comment As String
 
     Public Overrides Function ToString() As String
        If String.IsNullOrEmpty(Comment) Then
            Return String.Format("{0}={1}", IIf(Len(Title) > 0, Title, Name), Value)
        Else
            Return String.Format("{0}={1}; //{2}", IIf(Len(Title) > 0, Title, Name), Value, Comment)
        End If
    End Function
End Class

方程

一个方程代表一个生化反应，也表示底物（或系统组分）之间的相互作用。S-方程是一种化学速率方程，用于表示底物的浓度变化率。因此，方程类可以基本定义为两个属性

Public Class Equation
    <Xml.Serialization.XmlAttribute> Public Name As String
    <Xml.Serialization.XmlAttribute> Public Expression As String
End Class

它有一个计算方程表达式值的方法

Public Function Elapsed() As Boolean
    Call Me.sBuilder.Clear()
    sBuilder.Append(Expression)
    For Each e In Kernel.Vars 'Replace the name using the value
        Call sBuilder.Replace(e.Name, e.Value)
    Next
    Var.Value += (Microsoft.VisualBasic.Mathematical.Expression.Evaluate(sBuilder.ToString) * 0.1)
    Return True
End Function

声明

Var.Value += (Microsoft.VisualBasic.Mathematical.Expression.Evaluate(sBuilder.ToString)
* 0.1)

在此函数中，只进行迭代计算：使用方程表达式计算底物变化率，然后乘以经过时间以获得浓度变化值，最后使用变化值修改当前底物浓度。

扰动

生化系统的扰动是一种重要的研究方法，它可以让我们找出系统网络中的重要途径或反应。我们还可以通过比较模型中的计算结果或扰动结果与实验室实验结果，来检查某个途径或反应是否存在。

Public Class Disturb
    Public Enum Types
        Increase
        Decrease
        ChangeTo
    End Enum
    ''' <summary>
    ''' The name Id of the target.
    ''' (目标的名称)
    ''' </summary>
    ''' <remarks></remarks>
    <Xml.Serialization.XmlAttribute> Public Id As String
    ''' <summary>
    ''' The start time of this disturb.
    ''' (这个干扰动作的开始时间)
    ''' </summary>
    ''' <remarks></remarks>
    <Xml.Serialization.XmlAttribute> Public Start As Double
    ''' <summary>
    ''' The interval ticks between each kick.
    ''' (每次干扰动作执行的时间间隔)
    ''' </summary>
    ''' <remarks></remarks>
    <Xml.Serialization.XmlAttribute> Public Interval As Double
    ''' <summary>
    ''' The counts of the kicks.
    ''' (执行的次数)
    ''' </summary>
    ''' <remarks></remarks>
    <Xml.Serialization.XmlAttribute> Public Kicks As Integer
    <Xml.Serialization.XmlAttribute> Public DisturbType As Types
    <Xml.Serialization.XmlAttribute> Public Value As Double
End Class

在系统运行时，扰动对象对整个系统做了什么？扰动在特定时间直接修改底物浓度，因为底物浓度可以代表整个系统的状态。扰动有3种修改方法：

Public Shared ReadOnly Methods As System.Func(Of Double, Double, Double)() = {
    Function(Var As Double, Delta As Double) Var + Delta,
    Function(Var As Double, Delta As Double) Var - Delta,
    Function(Var As Double, Delta As Double) Delta}

    Public Enum Types
        Increase
        Decrease
        ChangeTo
    End Enum

扰动由Kicks对象管理，它有两个列表对象，一个用于未运行的扰动，另一个用于正在运行的扰动对象。

    ''' <summary> 
    ''' Standing by
    ''' </summary>
    ''' <remarks></remarks>
    <Xml.Serialization.XmlIgnore> Friend PendingKicks As List(Of Disturb)
    ''' <summary>
    ''' Active object. 
    ''' </summary>
    ''' <remarks></remarks>
    <Xml.Serialization.XmlIgnore> Friend RunningKicks As List(Of Disturb)

构建系统模块

脚本定义格式：所有声明语法均以“<关键字> 表达式”格式书写，注释行是指行首单词不在关键字列表中的行

关键字	信息	语法	示例
RXN	声明一个 S-方程	RXN <底物>=<表达式>	RXN X1=10X5-5X1^0.5
FUNC	声明一个函数	FUNC <名称> <表达式>	FUNC f x+y^2
INIT	设置底物的初始浓度值	INIT <底物>=<值>	INIT X1=2
CONST	声明一个常量	CONST <名称> <值>	CONST beta1 30
NAMED	设置底物的显示标题（可选）	NAMED <底物> <标题>	NAMED X1 ATP
STIMULATE	在系统中设置一个扰动	STIMULATE OBJECT <底物> START <开始时间> KICKS <踢次数> INTERVAL <踢间隔> VALUE <[类型]值>	STIMULATE OBJECT X3 START 45 KICKS 30 INTERVAL 7 VALUE 5
FINALTIME	设置系统运行时间	FINALTIME <值>	FINALTIME 10
TITLE, COMMENT	模型属性定义，可选	<关键字> <VALUE>	TITLE EXAMPLE SCRIPT

我在这里定义了一个Script类来将脚本文件编译成可用于内核的模型

这是一个脚本文件示例

/* This is a example script */
/* S-Equations */
RXN X1=10*X5*X3^-0.8-5*X1^0.5
RXN X2=5*X1^0.5-10*X2^0.5
RXN X3=2*X2^0.5-1.25*X3^0.5
RXN X4=8*X2^0.5-5*X4^0.5
// FUNC f 0.2*x+y^0.8
/* Substrate initial value */
// ATP comments ......
INIT X1=2
INIT X2=0.25
INIT X3=0.64
// BLA BLA BLA COMMENTS ......
INIT X4=0.64
INIT X5=0.5
/* Stimulates */
/* Kick  specific
times with add a value */
// Example comment 1
STIMULATE OBJECT X1 START 4 KICKS 1 INTERVAL 1 VALUE 0
// Example comment 2
/* STIMULATE OBJECT X2 START 5 KICKS -1 INTERVAL 2.5
VALUE --2 */
// NULL
/* STIMULATE OBJECT X3 START 45 KICKS 30 INTERVAL 7 VALUE
5 */
NAMED X1 ATP
NAMED X2 G6P
TITLE EXAMPLE SCRIPT
COMMENT NO COMMENT
FINALTIME 10

此脚本模拟了这样一个分支和反馈路径

设计模拟器内核

内核类很小，易于扩展，一个用于保存系统状态的数组，一个用于改变系统状态的表达式数组，一个用于执行扰动的模块，一个用于收集输出结果的模块，这四个组件几乎构成了一个完整的内核。

''' <summary>
''' The simulation system kernel.
''' </summary>
''' <remarks></remarks>
Public Class Kernel
    ''' <summary>
    ''' The current ticks since from the start of running.
    ''' (从运行开始后到当前的时间中所流逝的内核循环次数)
    ''' </summary>
    ''' <remarks></remarks>
    Friend RtTime As Double
    Friend Model As Model
    ''' <summary>
    ''' Data collecting
    ''' </summary>
    ''' <remarks></remarks>
    Dim DataAcquisition As New DataAcquisition
    ''' <summary>
    ''' Object that action the disturbing
    ''' </summary>
    ''' <remarks></remarks>
    Public Kicks As Kicks
    ''' <summary>
    ''' Store the system state.
    ''' </summary>
    ''' <remarks></remarks>
    Public Vars As Var()
    ''' <summary>
    ''' Alter the system state.
    ''' </summary>
    ''' <remarks></remarks>
    Public Channels As Equation()

内核有一个名为Run的方法，用于执行生化系统模拟：从零开始累积系统运行时间到最终时间，并在每个循环中进行迭代计算。是的，实际上整个系统状态的变化就是每个循环之间的迭代。

Public Sub Run()
    For Me.RtTime = 0 To Model.FinalTime Step 0.1
        Call Tick()
    Next
End Sub

''' <summary>
''' The kernel loop.(内核循环)
''' </summary>
''' <remarks></remarks>
Public Sub Tick()
    Dim Query As Generic.IEnumerable(Of Boolean) = From e As Equation In Channels Select e.Elapsed '
    Call DataAcquisition.Tick()
    Call Kicks.Tick()
    Query = Query.ToArray
End Sub

事实上，您还可以获得此内核的并行计算版本，只需简单修改LINQ查询语句即可：

Dim Query As Generic.IEnumerable(Of Boolean) = From e As Equation In Channels.AsParallel Select e.Elapsed '

但并行并不总是能提高效率，有时它可能会减慢系统运行速度，因为系统需要处理更多的事情。或者并行有时可能会导致计算数据不够理想。

输出仿真结果

内核获取一个数据采集类对象来收集模拟结果，并使用CSV文件存储数据。您可以在命名空间Microsoft.VisualBasic.DataVisualization.Csv.File中找到CSV包装器类。

Public Class DataAcquisition 
    Dim _dataPackage As New Microsoft.VisualBasic.DataVisualization.Csv.File
    Dim Kernel As Kernel

    Public Sub Tick()
        Dim Row As New Microsoft.VisualBasic.DataVisualization.Csv.File.Row
        Call Row.Add(Kernel.RtTime)
        For Each Var As Var In Kernel.Vars
            Row.Add(Var.Value)
        Next
        _dataPackage.AppendLine(Row)
    End Sub

    Public Sub Save(Path As String)
        Call _dataPackage.Save(Path)
    End Sub

    Public Sub [Set](Kernel As Kernel)
        Dim Row As New Microsoft.VisualBasic.DataVisualization.Csv.File.Row
        Me.Kernel = Kernel
        Call Row.Add("Elapsed Time")
        For Each Var As Var In Kernel.Vars
            If Not String.IsNullOrEmpty(Var.Title) Then
                Call Row.Add(String.Format("""{0}""", Var.Title))
            Else
                Call Row.Add(Var.Name)
            End If
        Next
        Call _dataPackage.AppendLine(Row)
    End Sub

    Public Shared Function [Get](e As DataAcquisition) As Microsoft.VisualBasic.DataVisualization.Csv.File
        Return e._dataPackage
    End Function
End Class

测试

这里我使用一个示例系统来测试程序：Hull et al 基准测试 / 非刚性系统。

系统定义如下

RXN z1=2*(z1-z1*z2)
RXN z2=-1*(z2-z1*z2)
INIT z1=2
INIT z2=3
FINALTIME 20

将此模型放入文本文件，然后使用代码运行此系统

Sub Main()
    Dim CSV = Kernel.Run(Script.Compile("./Hull.txt"))
    CSV.Save("/home/xieguigang/Documents/Hull.csv")
End Sub

然后我们打开从计算中获得的CSV数据文件，在WPS Office或OpenOffice中打开，使用数据绘制散点图。

从模拟结果中我们可以看到，两种底物z1和z2的浓度随时间周期性变化，并且具有良好的周期性。

关注点

s-系统利用线性方程组模拟非线性复杂系统。虽然复杂系统被我们简化了，但我们仍然可以通过模拟发现一些有趣的属性：复杂系统对每个组件的初始值敏感。

我们可以改变hull系统中组件z1的初始值，然后运行修改后的模型，看看我们得到了什么

Z1=1

Z1=0

Z1=-1

从z1初始值的变化中我们可以看到，系统开始失去在周期时间内改变组件值的能力。最有趣的是，当z1=1时，系统在第一时间具有周期性属性，但随后情况变得疯狂：持续的系统计算误差累积最终改变了整个系统行为。我们也知道这种现象的著名名称：蝴蝶效应，累积的系统误差最终会改变整个系统行为。

生化系统也做同样的事情吗？

尽管生命系统也是一个复杂的系统，但如果干扰不那么致命，生命系统具有从干扰中恢复的能力（我们对系统初始值的修改是对正常状态的干扰）。系统生物学指出，生命系统具有以下特性：涌现性、鲁棒性、复杂性和模块化。因此，从这些角度来看，我们可以得出结论，如果进行计算，生物生命系统会产生大量的系统误差，并且由于其鲁棒性，生物系统对生化参数的精确值不敏感，它可以从误差中自我调整，但数学方程不能。为什么生物系统不那么敏感？因为它已经找到了从系统误差中恢复或避免系统崩溃的方法，即调控网络，它确实非常复杂。调控网络是一个真正的复杂系统，而数学方程是一个虚假的系统，这就是原因。

从这个角度来看，这意味着方程不足以模拟一个复杂系统。虽然数学方程可以给我们数值结果，但这个数值结果只是一个简化复杂系统的计算结果。

因此，基于E.O.Voit的基础工作以及S-系统的概念和工作机制，我们创建了一个PLAS系统的升级版本，称之为Object-S。该系统不仅包含数学方程，还整合了生物大分子对象、相互作用关系和细胞事件。我将这个新的系统生物学建模项目称为基因时钟计划（GCI），正如微软也有一个名为Microsoft Biology Initiative (MBI) 或 Microsoft Biology Foundation (MBF) 的.NET生物信息学研究项目。

该程序在Ubuntu 13.04 (mono 2.1) 上开发，并成功在Ubuntu LINUX 和 Windows 8 企业版上调试和测试。