摘要：MATLAB是美国Mathwork公司于20世纪80年代中期推出的适用于科学和工程计算的数学软件。它具有优秀的数值计算能力、卓越的数据可视化功能和丰富的库函数资源等优点，使其在数学软件中脱颖而出。本文以生化反应的动态变化为例，用MATLAB模拟计算生化反应的动态过程。
关键词：生化反应；模拟；动态过程；MATLAB
中图分类号：O563.9 文献标示：A
引言
生化反应数学模型一般比较复杂，求解不太容易，浓度变化过程不太好直观表示出来。应用数学理论和数值方法定性或定量地预测或模拟各物种在生化反应过程中的浓度随时间的变化关系是反应动力学研究的主要内容之一。确定浓度随时间的变化关系也就是求解各反应组分与时间的微分关系问题。计算机技术的发展为模拟计算生化反应动态过程起到了巨大的推动作用。
MATLAB以数组和矩阵为基础，最初主要用于数值计算和自动控制领域，现已逐步拓展为数据可视化、数字信号处理、图像处理、数理统计和经济等领域。有文献[1-3]介绍了MATLAB在化工模拟计算中的简单应用。本文着重介绍MATLAB在生化反应模拟计算中的应用，突出表现MATLAB丰富的函数库和数据可视化功能。
求解非刚性常微分方程组的函数
一阶常微分方程（组）初值问题包括刚性和非刚性问题两大类，但化工领域中多数情况属于非刚性ODE问题。MATLAB提供了的非刚性方程积分函数有以下三个[4]：
ode45——适用于精度较高的四五阶Runge-Kutta法；
ode23——适用于精度较低的二三阶Runge-Kutta法；
ode113——适用于可变阶dams-Bashforth-Moulton法。
计算机性能很好，运行速度极快，一般尽量适用精度较高的算法。非刚性常微分方程组常用ode45函数求解。
Ode45函数的用法：
[t, y] = ode45(@ODEfun,tspan,y0)
[t,y]=ode45(@ODEfun, tspan,y0,options,p1,p2…)
其中，ODEfun为自定义的函数名，内容是微分方程组的矩阵表带；tspan为求值范围，可以是积分限，也可以是一些离散点的向量；y0为状态变量的初始条件向量（一般是列向量）；options为积分参数选项（可选），options=[]时表示取默认值；p1, p2…为直接传递给自定义函数ODEfun的已知参数；t, y为输出变量，在生化反应中一般为时间和浓度。
MATLAB用于生化反应动态、
静态模拟计算的实现
● 用ode45函数模拟厌氧间歇发酵动态过程
使用Saccharomyces cerebisiae（一种干酵母）厌氧间歇发酵葡萄糖生成乙醇和副产物二氧化碳。物料包含10%的葡萄糖。接种体是0.0005（质量比）。在14%乙醇时，由于产品抑制使细胞停止生长。假设kd=0但当底物完全消耗时，开始消耗细胞。模拟计算从发酵开始各组分的浓度动态变化[5]。已知初始浓度为：t=0,X0=0.0005,S0=0.04,p0=0。已知参数为：umax=0.5h-1, Ks=0.001, m=1, Ms=0.008h-1, ^Yx/s=2, Yx/s=1。
数学模型：
dX/dt=Rs (1)
dS/dt=Rs (2)
dp/dt=Rp (3)
Rx=umaxX[S/(Ks+S)](1-p/pmax)m (4)
Rs=-Rx/Yx/s-MsX (5)
Rp=Rx/Yx/s-Rx/^Yx/s+MsX (6)
方程（1）、（2）、（3）分别为细胞、底物和产物的质量平衡方程;(4)、(5)、(6) 分别为抑制时细胞生长速率方程、底物消耗的速率方程和产物生成的速率方程。
根据上面的数学模型，建立求解微分方程组的函数：
global umax KS pmax m YX2S YX2Shat MS
%已知参数和初始浓度
umax=0.5; pmax=0.109; KS=0.001; m=1; YX2Shat=2; YX2S=1; MS=0.008; X0=0.0005; S0=0.04; p0=0;
%ode45( )求解微分方程组
tspan=[0 15]; y0=[X0 S0 p0];
[t y]=ode45(@ModelE,tspan,y0)
%绘制浓度动态变化图
plot(t,y, 'k'), xlabel('反应时间/h'), ylabel('质量分数');gtext('活细胞'),gtext('葡萄糖'), gtext('乙醇');
%模型方程
function dydt=ModelE(t,y,k)
global umax KS pmax m YX2S YX2Shat MS
X=y(1); S=y(2); p=y(3);
RX=umax×X×S/(KS+S)×(1-p/pmax)^m;
RS=-RX/YX2S-MS×X;
Rp=RX/YX2S-RX/YX2Shat+MS×X;
dydt=[RX;RS;Rp];
执行程序后，得到动态变化过程如图1所示。
3.2 用ode45和fsolve函数模拟计算生化反应器[6]
在一个连续搅拌槽式发酵罐中进行生化反应，小x1为生物量的浓度，x2为酶解物（底物）的浓度，xf1、xf2分别为进料中所含的生物量和酶解物的浓度，F为料液的体积流量，V为工作体积。已知数据为：D=0.3，umax=0.53，Y=0.4，km=0.12，xf1=0，xf2=0.4，k1=0.4545。对生物量和底物进行稳态、动态模拟[6]。
数学模型：
动态时， dx1/dt=(u-D)x1 (7)
dx2/dt=D(xf2-x2)-ux1/Y (8)
稳态时， (u-D)x1=0 (9)
D(xf2-x2)-ux1/Y=0 (10)
式中 D=F/V，稀释速度，s-1； u=umaxx2/(km+x2)或u=umaxx2/(km+x2+k1x22)，生长速率常数；km和umax为动力学参数；Y为收率；x1和x2为浓度；xf1和xf2表示进料时浓度。初始条件为，t=0，x1=1，x2=1。
根据数学模型方程（7）～（10），编写求解微分方程组的函数：
global D umax Y km xf1 xf2 k1 ;
%已知参数和变量
D=0.3; umax=0.53; Y=0.4; km=0.12; xf1=0; xf2=4.0; k1=0.4545;
%ode45( )求解微分方程
tspan=[0 25]; y0=[1 1];
[t y]=ode45(@ModelDyn,tspan,y0);
%绘制动态变化图
plot(t,y);xlabel('反应时间/h');ylabel('浓度x');gtext('生物量'),gtext('底物');
%fsolve( )求解稳态方程组
x0=[1 1]'; x=fsolve(@ModelSta,x0);
%动态时模型方程
function dydt= ModelDyn(t,y,k)
global D umax Y km xf1 xf2 k1 ;
x(1)=y(1);x(2)=y(2);
u=umax×x(2)/(km+x(2));
%或u=umax×x(2)/(km+x(2)+k1×x(2)^2)
Rx1=(u-D)×x(1);
Rx2=D×(xf2-x(2))-u×x(1)/Y;
dydt=[Rx1;Rx2];
%稳态时模型方程
function f=ModelSta(x)
global D umax Y km xf1 xf2 k1 ;
u=umax×x(2)/(km+x(2));
%或u=umax×(2)/(km+x(2)+k1×x(2)^2)
f(1)= (u-D)×x(1);
f(2)= D×(xf2-x(2))-u×x(1)/Y;
执行程序结果为：当u=umaxx2/(km+x2)时，生物量和底物浓度的稳态值分别为0.0和4.0，浓度动态变化如图2所示。当u=umaxx2/(km+x2+k1x22)时，生物量和底物浓度的稳态值分别为0.9951和1.5122，浓度动态变化如图3所示。
结语
从上面的模拟实现可以看出，利用MATLAB丰富的库函数资源可以让编程人员从繁琐的程序代码中解放出来，并且提供的数据可视化功能，所以特别适合于工程人员的计算模拟。MATLAB是新一代计算机语言，程序编写方便、简洁且极易调试，人机交互性强。因此，在化学化工领域模拟计算中，MATLAB也会发挥其特长，必将成为化工专业人员在模拟计算中的有力工具。
关键词：生化反应；模拟；动态过程；MATLAB
中图分类号：O563.9 文献标示：A
引言
生化反应数学模型一般比较复杂，求解不太容易，浓度变化过程不太好直观表示出来。应用数学理论和数值方法定性或定量地预测或模拟各物种在生化反应过程中的浓度随时间的变化关系是反应动力学研究的主要内容之一。确定浓度随时间的变化关系也就是求解各反应组分与时间的微分关系问题。计算机技术的发展为模拟计算生化反应动态过程起到了巨大的推动作用。
MATLAB以数组和矩阵为基础，最初主要用于数值计算和自动控制领域，现已逐步拓展为数据可视化、数字信号处理、图像处理、数理统计和经济等领域。有文献[1-3]介绍了MATLAB在化工模拟计算中的简单应用。本文着重介绍MATLAB在生化反应模拟计算中的应用，突出表现MATLAB丰富的函数库和数据可视化功能。
求解非刚性常微分方程组的函数
一阶常微分方程（组）初值问题包括刚性和非刚性问题两大类，但化工领域中多数情况属于非刚性ODE问题。MATLAB提供了的非刚性方程积分函数有以下三个[4]：
ode45——适用于精度较高的四五阶Runge-Kutta法；
ode23——适用于精度较低的二三阶Runge-Kutta法；
ode113——适用于可变阶dams-Bashforth-Moulton法。
计算机性能很好，运行速度极快，一般尽量适用精度较高的算法。非刚性常微分方程组常用ode45函数求解。
Ode45函数的用法：
[t, y] = ode45(@ODEfun,tspan,y0)
[t,y]=ode45(@ODEfun, tspan,y0,options,p1,p2…)
其中，ODEfun为自定义的函数名，内容是微分方程组的矩阵表带；tspan为求值范围，可以是积分限，也可以是一些离散点的向量；y0为状态变量的初始条件向量（一般是列向量）；options为积分参数选项（可选），options=[]时表示取默认值；p1, p2…为直接传递给自定义函数ODEfun的已知参数；t, y为输出变量，在生化反应中一般为时间和浓度。
MATLAB用于生化反应动态、
静态模拟计算的实现
● 用ode45函数模拟厌氧间歇发酵动态过程
使用Saccharomyces cerebisiae（一种干酵母）厌氧间歇发酵葡萄糖生成乙醇和副产物二氧化碳。物料包含10%的葡萄糖。接种体是0.0005（质量比）。在14%乙醇时，由于产品抑制使细胞停止生长。假设kd=0但当底物完全消耗时，开始消耗细胞。模拟计算从发酵开始各组分的浓度动态变化[5]。已知初始浓度为：t=0,X0=0.0005,S0=0.04,p0=0。已知参数为：umax=0.5h-1, Ks=0.001, m=1, Ms=0.008h-1, ^Yx/s=2, Yx/s=1。
数学模型：
dX/dt=Rs (1)
dS/dt=Rs (2)
dp/dt=Rp (3)
Rx=umaxX[S/(Ks+S)](1-p/pmax)m (4)
Rs=-Rx/Yx/s-MsX (5)
Rp=Rx/Yx/s-Rx/^Yx/s+MsX (6)
方程（1）、（2）、（3）分别为细胞、底物和产物的质量平衡方程;(4)、(5)、(6) 分别为抑制时细胞生长速率方程、底物消耗的速率方程和产物生成的速率方程。
根据上面的数学模型，建立求解微分方程组的函数：
global umax KS pmax m YX2S YX2Shat MS
%已知参数和初始浓度
umax=0.5; pmax=0.109; KS=0.001; m=1; YX2Shat=2; YX2S=1; MS=0.008; X0=0.0005; S0=0.04; p0=0;
%ode45( )求解微分方程组
tspan=[0 15]; y0=[X0 S0 p0];
[t y]=ode45(@ModelE,tspan,y0)
%绘制浓度动态变化图
plot(t,y, 'k'), xlabel('反应时间/h'), ylabel('质量分数');gtext('活细胞'),gtext('葡萄糖'), gtext('乙醇');
%模型方程
function dydt=ModelE(t,y,k)
global umax KS pmax m YX2S YX2Shat MS
X=y(1); S=y(2); p=y(3);
RX=umax×X×S/(KS+S)×(1-p/pmax)^m;
RS=-RX/YX2S-MS×X;
Rp=RX/YX2S-RX/YX2Shat+MS×X;
dydt=[RX;RS;Rp];
执行程序后，得到动态变化过程如图1所示。
3.2 用ode45和fsolve函数模拟计算生化反应器[6]
在一个连续搅拌槽式发酵罐中进行生化反应，小x1为生物量的浓度，x2为酶解物（底物）的浓度，xf1、xf2分别为进料中所含的生物量和酶解物的浓度，F为料液的体积流量，V为工作体积。已知数据为：D=0.3，umax=0.53，Y=0.4，km=0.12，xf1=0，xf2=0.4，k1=0.4545。对生物量和底物进行稳态、动态模拟[6]。
数学模型：
动态时， dx1/dt=(u-D)x1 (7)
dx2/dt=D(xf2-x2)-ux1/Y (8)
稳态时， (u-D)x1=0 (9)
D(xf2-x2)-ux1/Y=0 (10)
式中 D=F/V，稀释速度，s-1； u=umaxx2/(km+x2)或u=umaxx2/(km+x2+k1x22)，生长速率常数；km和umax为动力学参数；Y为收率；x1和x2为浓度；xf1和xf2表示进料时浓度。初始条件为，t=0，x1=1，x2=1。
根据数学模型方程（7）～（10），编写求解微分方程组的函数：
global D umax Y km xf1 xf2 k1 ;
%已知参数和变量
D=0.3; umax=0.53; Y=0.4; km=0.12; xf1=0; xf2=4.0; k1=0.4545;
%ode45( )求解微分方程
tspan=[0 25]; y0=[1 1];
[t y]=ode45(@ModelDyn,tspan,y0);
%绘制动态变化图
plot(t,y);xlabel('反应时间/h');ylabel('浓度x');gtext('生物量'),gtext('底物');
%fsolve( )求解稳态方程组
x0=[1 1]'; x=fsolve(@ModelSta,x0);
%动态时模型方程
function dydt= ModelDyn(t,y,k)
global D umax Y km xf1 xf2 k1 ;
x(1)=y(1);x(2)=y(2);
u=umax×x(2)/(km+x(2));
%或u=umax×x(2)/(km+x(2)+k1×x(2)^2)
Rx1=(u-D)×x(1);
Rx2=D×(xf2-x(2))-u×x(1)/Y;
dydt=[Rx1;Rx2];
%稳态时模型方程
function f=ModelSta(x)
global D umax Y km xf1 xf2 k1 ;
u=umax×x(2)/(km+x(2));
%或u=umax×(2)/(km+x(2)+k1×x(2)^2)
f(1)= (u-D)×x(1);
f(2)= D×(xf2-x(2))-u×x(1)/Y;
执行程序结果为：当u=umaxx2/(km+x2)时，生物量和底物浓度的稳态值分别为0.0和4.0，浓度动态变化如图2所示。当u=umaxx2/(km+x2+k1x22)时，生物量和底物浓度的稳态值分别为0.9951和1.5122，浓度动态变化如图3所示。
结语
从上面的模拟实现可以看出，利用MATLAB丰富的库函数资源可以让编程人员从繁琐的程序代码中解放出来，并且提供的数据可视化功能，所以特别适合于工程人员的计算模拟。MATLAB是新一代计算机语言，程序编写方便、简洁且极易调试，人机交互性强。因此，在化学化工领域模拟计算中，MATLAB也会发挥其特长，必将成为化工专业人员在模拟计算中的有力工具。
