Loading

MODELING

1.Summary

Recent research proves, NO-cGMP pathway has important significance to angiocardiopathy's treatment. Target related to the pathway has recently been widely discussed in drug research and development. As members of Modeling team of iGEM in South China University of Technology, we hope to explore NO-cGMP pathway on the aspects of mathematical theory and model.

We analysed the pathway's primary link's situation,including using ODE to analyse the process of NO activating sGC and Michealis Equation to analyse cGMP's synthesis and degradation.

According to the differential equation derived from the theory, we got the relationship between the yield of sGC in different states,the yield of cGMP, the production and degradation rate of cGMP and time,and the relationship between the yield of bioactive sGC, cGMP and NO by using Matlab programing. In addition, we studied the effect of different feedback inhibition levels of cGMP on the amount of the above, and the relationship was used as the feasible criteria for the study of this pathway.

Due to the limitations of our students' own level, the process of producing PKG in the cGMP can only be based on the experience of the usual mathematical modeling, using the normalized method and Hill Equation for a little analysis.

In order to study the relationship between the amount of cGMP and sGC, we also use Matlab to draw the image of the relationship by using the fitting function, so as to have a more intuitive understanding of the process.

Although scientists are not thorough enough about the study of Lentivirus, the iGEM modeling group still wants to use the existing model to establish a mathematical model to describe the process of gene transfection and combination in cells.

2.Theory

2.1 NO activates sGC

NO activates sGC from substrate state D1 to activated state D2. sGC is a kind of heterodimer whose center is composed of two subunits, which are called the α and β subunits of ferroheme. NO stimulates sGC by combining sGC’s ferroheme(300-700 times). Recent research shows the combining mechanism is involved with two steps: During the first step, NO combines substrate sGC to generate the ferrous nitrosyl heme complexes D3 which then decays to a fully activated state mixture D2.

The first step of NO combined with sGC can be expressed as a reversible biochemical reaction process:

D1 and D3 represent the initial and intermediate states of sGC. Accordingly, k and k-1 are rate constants. The second step has two parallel paths which are depend or independ on NO. The two are all involved in the transition from D3 to D2. The process of independent NO was considered to be a natural decay from the state D3 to D2, with the His105 residue separated from the heme. The process can be expressed as a first order and irreversible kinetic equation(rate constant is k2):

The path which depends on NO can be described as enzymatic action based on the traditional Michaelis Menten kinetics. Since the reversible and dissociation reactions of NO and sGC heme were rapid (association rate constant is 7.1±2*10⁸ M^-1s^-1 in 4℃), the damage of His105 residue becomes the rate limiting step in the process (the observed value of rate constant of the process which depends on NO in reference documentation is 2.4*10⁵M^-1s^-1 in 4℃). Therefore, we assume that the NO dependent pathway can be equivalent to a single step irreversible process with a rate constant k3:

Another side of NO-sGC reaction’s is NO disaggregates from sGC, which means D2 is inactivated、sGC return to the initial state. We can express the process as a first order and irreversible kinetic equation:

Rate constant k4 can be regulated by multiple factors, including NO scavenger, substrate concentration and so on. Under the stimulation of cGMP, the yield of NO presents biphase. Along with the initial rapid growth to the peak，following is a recovery to a lower level, with the increase in degradation rate makes the concentration of cGMP reaches to a steady state. Bellamy suggested that the unknown sGC regulatory factors may exist and that the sensitivity of sGC to NO in the late state of cGMP reaction is decreased, and these regulatory factors are active in response to NO's exogenous and endogenous. Now, we assume that the local feedback loop of the final product's cGMP level is a regulator of sGC desensitization. In the hypothesis that cGMP may not be a direct inhibitor of sGC activity, it may be the second messenger that generates negative feedback. We assume that NO separated from sGC and D2 becoming D1 are cGMP mediated processes. In our model, the rate constant of D2 returning to D1 can be expressed as:

k₄=K₄[cGMP]^m

K4 is constant. Here, we assume that k4 is mth order of cGMP concentration, where m is the feedback strength of cGMP.

Summarizing the analysis of sGC kinetics above, we using differential equation to describe the whole process:

Let D1=[D1]/[D0], D3=[D3]/[D0], D2=[D2]/[D0]. Meanwhile, we assume that the total concentration of sGC is constant, that is to say D1+D3+D2=1.

We suppose that the NO is immediately consumed when NO is separated from the D2 or leaves the sGC, and generate other substances that make NO no longer free but produce other biochemical components. Thus the result of NO-sGC interaction is that NO have a net loss. The total consumption rate of NO including NO binding sGC can be approximated as a first order kinetic equation:

Here, Lno is the total amount of NO in cells and in the outside, which represens the available amount NO in physiological or experimental conditions, t is time; Kno express the total consumption rate constant of NO, which reflects the bioactivity of NO.

2.2 cGMP's synthesis and degratation

Since the synthesis of cGMP is after sGC(D2)'s forming which has catalytic activity, we assume the yield of cGMP generated from GTP mainly depends on the degree of D2’s activity. cGMP's life is determined by PDE5A's multiple subtypes. PDE5A can be hydrolyzed to GMP. The process of production and degradation can be modeled by Michealis equation:

Supposing that GTP is abundant, the yield of cGMP mainly depend on the concentration of sGC. Studies suggested PDE5A with different substrate affinity、subcellular's position and kinetic properties can catalysis cGMP's hydrolysis. cGMP's balance equation is:

d[cGMP]/dt = V_sGC,maxD2- V_PDE,max[cGMP]/(K_m,PDE+[cGMP])

Inside,V_sGC,max represents cGMP's maximum production rate when D2=1, which means all avaliable enzyme molecules participated in catalysis.V_PDE,max is cGMP's maximum hydrolyzation rate.V_sGC,maxis directly proportional to sGC's total concentration[D0], and D0 is constant in our model. Consequently, cGMP's production rate and degradation rate are:

V_p=V_sGC,maxD2

and

Michaelis constant K_m,PDE is measured in experiment by Turko et al. Recent studies show, during the process of hydrolyzation the basic level of PDE5A's bioactivation is relatively low. However, after cGMP combining NH2-end structure regulation by itself, this bioactivation is greatly improved(up to 9-11 times). Through cGMP, PDE5A's activated process can be discribed as PDE_inact + cGMP↔ PDE_act. Therefore, active PDE5A's concentration depends on the level of cGMP which because V_PDE,max is directly proportional to active PDE5A's total concentration in Michaelis kinetics. Meanwhile, we suggest that active PDE5A's concentration is nearly proportional to cGMP's concentration. It's a reasonable approximation, if active PDE5A's level is relatively high, cGMP's level is relatively high accordingly. Thus, we consider cGMP's maximum hydrolyzation V_PDE,max is directly proportional to [cGMP]:

V_PDE,max=k_PDE[cGMP]

table 1

Model Parameters
Parameter	Value
k₁	2.0nM^-1s^-1
k_-1	100s^-1
k₂	0.1s^-1
k₃	0.003nM^-1s^-1
K_m,PDE	2.0μM
K_no	0.01s^-1

table 2

Model Parameters
m	V_sGC,max,μM	K_PDE,s^-1	K₄
0	1.26	0.0695	0.4
1	1.09	0.032	0.98
2	0.8520	0.0195	0.011

Above is model's parameters from past published documents: these parameters' values are all determined under the condition of NO's amount of substance is 220nH, at this time，D1=1，D2=D3=0，cGMP’s initial concentration is 0，V_sGC,max, K_PDE, K₄is constant.

2.3 Innovation analysis of cGMP combined with PKG

In our modeling, we assume that, compared with the synthetic rate of cGMP, PKG was activated by cGMP and the target cell phosphorylationis is relatively a fast event. The regulation of cGMP can be evaluated by the normalization function of the concentration of cGMP, A ([cGMP]): when the available amount of cGMP reached the minimum (i.e.: [cGMP]→0, namely: no impact), A→0 ; when the available amount of cGMP reached the maximum (i.e.: the greatest impact), A→1. Therefore, we consider cGMP as a molecular switch of adjusting start and termination of the reaction. We use the Hill Equation to approximate the function:

K_m and the Hill coefficient(n_H) relatively represent the value when variation tendency of model changes and the step degree of switch. This is only an approximate empirical formula, we can get the precise equation by experiment data. Considering that the reaction rate of cGMP combined with PKG is fast, we assume that the regulation process is reasonable. Therefore, any value of target number Q which is generated by the changes of the amount of cGMP can be generally described by the following algebraic equation:

Q=Q_basal+A_cGMP(K_m,n_H)Q_change

Q_basal represents the basic situation of cGMP, Q_changeepresents the biggest change of which is induced by cGMP. Q can represent the rate constant of enzymatic reaction, conductivity of ion channels and so on.

But due to the restrictions of our students’ own level, we can’t make mathematical reasoning for the complex process deeply, so we can only write the main idea as above. As to describe this process by mathematical model language completely and to draw a conclusion which is consistent with the actual situation, we can complete it only to wait until after we improve ourselves.

3.Results Analysis (sGC desensitization: cGMP local feedback inhibition)

figure.1

figure.2

figure.3

By Figure.1 ,we see that in the linear feedback conditions of cGMP(i.e.: M = 1), different state of sGC changes with the time of change: D1 plummets in a brief period when the reaction started, and then reached a nadir in 8 seconds or so. After a short rising stage it tends to be stable finally; D3, the intermediate product of reactions of NO and sGC experienced a period of increase in the first time and then went to the attenuation and finally reached the stability; D2 (i.e.: high-active sGC) reachs the steady state as the same after a brief surge and small amplitude fluctuations. Figure.2 shows the production rate of cGMP experienced process of rising and falling and the degradation rate of cGMP changes from fast to slow and how they tends to be uniform eventually. Figure.3 shows the phenomenon that the production of cGMP experienced process of rising and eventually tends to be stabilized.

Therefore, the autoregulation of cells has the function of making every sector of NO-cGMP pathway tend to be stabilized.

figure.4

figure.5

figure.6

figure.7

In order to discuss how the feedback inhibition of cGMP influences the desensitization of sGC and the degree, we use Matlab to make four images, namely, Figure.4 to Figure.7. Figure.4 and Figure.5 are images of no feedback inhibition (i.e.: m=0). The changes of D1, D2, D3, Vp, Vd with time are very simple, either increase or decrease to stabilization, and there is no volume fluctuate with time during the period. However, when the value of the feedback degree of m is 2, great changes have occurred. Fluctuation degree of each volume is more intense compared with all sectors when m=1. But without any exception, all the volumes tend to be stable, which shows that our model is consistent with the actual law.

figure.8

figure.9

Figure.8 and Figure.9 respectively show the relationship between NO and cGMP along with sGC in the case of cGMP linear feedback inhibition (m=1). We find that in feedback inhibition of cGMP the volume of cGMP is always less than that with no feedback inhibition. Meanwhile, apart from appearing sharply rising in the beginning and is higher than that with no feedback inhibition, as time goes by, the amount of sGC of feedback inhibition gradually reduces and is lower than that with no feedback inhibition and hold it on. In summary, the function of feedback inhibition of cGMP changes the relationship between cGMP-NO and sGC-NO into a lower state which is consistent with many research results.

4.Dosage effect of sGC on cGMP

cGMP(nmol/L)	185.5961	189.632	221.0417	287.9465	306.8128
sGC(U/L)	11.78422893	12.51607142	13.93817204	17.98381967	20.11388292

figure.10

In order to study the effect of sGC on cGMP, we need to find out the curve of the experimental group. Since the coordinates of the starting point in the coordinate system is (180,11) and the growth rate of cGMP with the increase of sGC concentration is more and more fast, we might as well set up the relationship between y and x as:

y=a+be^kx

By using ‘lsqcurvefit’ function in Matlab, we got the value of a=8.3573, b=0.6748, k=0.0093.

figure.11

From the image above, the bigger the concentration of sGC is, the bigger the concentration of cGMP, and the growth rate of cGMP is increasing. However, the concentration of cGMP can't increase without limit in the cell, but it will reach to a certain level of stability. This data is only from the experimental group's exploratory experiments which only reflects the change of cGMP in a certain range with the sGC. In order to draw a more accurate conclusion, more sophisticated experiments should be done.

5.Transfection and gene combinations of Lentivirus in HEK293 cells

iGEM experimental group used Lentivirus in the experiment. Lentivirus is gene therapy vectors which are developed on the basis of HIV-1,but so far, scientists have not researched the details of the HIV infection, integration, replication and so on clearly, it also means that lentivirus mechanism is not very clear. But the macro and micro mechanism of the exploration of Lentivirus is indeed attracting lots of people interested involved in the study. Here, the iGEM modeling group wants to use some of the existing models to establish the mathematical model to describe the process of the transfection and gene combinations of Lentivirus in the cells.

5.1 Macro mechanism

In order to explore whether Lentivirus has a time delay like HIV in the erosion of the cell, we first make a hypothesis. We assume there is a time delay between the initial Lentivirus invasion susceptible cells and the subsequent infection of the virus. Because the Lentivirus can not produce new virus particles, so when we analyze the whole process,we should avoid the generation of new virus particles. Therefore, according to the research results of Perelson、Bonhoeffer and Xu Rui, we have integrated to establish the model of Lentivirus’s transfection with time delay:

Among them, healthy cells (also known as susceptible cells) x(t) produce with constant input rate,d represents the death rate, represents saturated infection rate, healthy host cells are exposed to the virus particles after infected with bilinear infection rate of βxv into infected cells y(t);α>0; Infected cells’ death rate is a; The removal rate and initial amount of dissociative virus is μ and k; τ represents the time of virus’ invasion into susceptible cells to susceptible cells infected until death. Assuming m represents the production rate of infected cells(not dead), e^-mτdenotes the survival probability of the infected cells from time t-τ to time t (1/m represents the average survival time of cells infected but not with infectivity).

figure.12

In the case of other constant unchanged, we take different τ (τ=15,5,2.5,1,0.1) and use dde in MATLAB to solve the delay differential equations. From the image above, there are no changes of images with the changes of τ, but always overlap. Thus, we conclude that based on the establishment of Delay Lentiviral Transfection Model, Lentivirus does not have a time delay, which is different with HIV. As to whether Lentivirus completely has no time delay or it has other special properties, further studies need to be done in the future.

5.2 Micro mechanism

There is no conclusion on the study of Lentivirus, but because Lentivirus is gene therapy vectors which are developed on the basis of HIV-1, so we can simulate the process of cells infected by Lentivirus according to the process of HIV-1 infecting cells and some of the characteristics of Lentivirus. As Lentivirus belongs to RNA virus and the RNA virus’ coat protein contains the reverse transcriptase, the reverse transcriptase of Lentivirus after the invasion of the cells is brought into the cell by the virus itself. So there are differential equations as follows:

Reverse transcriptase:

The synthesis of cDNA:

The synthesis of double stranded DNA:

The synthesis of provirus DNA:

DNA negative strand transcription:

DNA positive strand replication:

Five kinds of nucleotides:

Among them, x1 to x5 represent positive strand, negative strand, cDNA, double stranded DNA and provirus DNA relatively, x6 to x10 represent five nucleotides, x11 represents the reverse transcriptase, bi(i=6,……,10)，ci(i=1,2)，di(i=1,……,5) and mi(i=1,2,3,4,11) are all constant. v1 and v2 are the maximum reaction rate for the two reactions; δ 1 represents the initial amount of reverse transcriptase in Lentivirus; δi(i=6,……,10) represent the five kinds of the initial amount of nucleotides and they are all constant.

figure.13

From the image above, the quantity of the intermediate products produced by Lentivirus in the process of making the host cell death increased rapidly and reached stable, which, in quantity, can be directly explained why Lentivirus can not produce the structural proteins as HIV virus and achieve the assembly of the virus particles, so that the amount of the intermediate production continued to grow.

6.Discussion

In the process of modeling, there exists a lot of problems. First, considering the production of cGMP can inhibit the reaction system. therefore, we assume that NO seperated from sGC and D2 transforming into D1 are the cGMP mediated process, so we use the feedback factor m to show the degree of the feedback. The question here is whether the degree of feedback will be in accordance with the order we have chosen. Second, there is no conclusion about the research of Lentivirus and the model is only based on the HIV research results of scientific workers. Therefore, in this sense, the model of Lentivirus’ transfection and gene combination in HEK293 cells still exist defects.

Despite the above limitations, our model can still an extensible framework for the study of the important NO-cGMP signaling mechanisms. Considering the big variety of cells across different tissues and species, NO-cGMP mechanisms are bound to vary and have designs that are specific to cell types. Thus,the model can serve as a general modeling framework for the study of NO-cGMP effects on the HEK293 as a complex cell signaling system. This framework can be tailored to address particular situations and can also be expanded to incorporate newly observed mechanisms.

In addition, based on scientific workers’ experimental and modeling research of HIV, we have established the model of simulating Lentivirus’ invading cells. Therefore, this model can be used as a reference for other iGEM teams to study Lentivirus.

Reference:

[1] A First Course in Mathematical Modeling(Third Edition) Frank R.Giordano, Maurice D.Weir, William P.Fox

[2] Mathematical Models Qiyuan Jiang, Jinxing Xie, Jun Ye (2010)

[3] Mo E, Amin H, Bianco IH, and Garthwaite J. Kinetics of a cellular nitric oxide/cGMP/phosphodiesterase-5 pathway. J Biol Chem 279: 26149–26158, 2004.

[4] Lucas KA, Pitari GM, Kazerounian S, Ruiz-Stewart I, Park J, Schulz S, Chepenik KP, and Waldman SA. Guanylyl cyclases and signaling by cyclic GMP. Pharmacol Rev 52: 375–414, 2000.

[5] Friebe A and Koesling D. Regulation of nitric oxide-sensitive guanylyl cyclase. Circ Res 93: 96–105, 2003.

[6] Global stabitility of infectious diseases models with saturation infection and humoral immunity Yonglu Tang,2012.

Appendix:

Source code

1.

1)When m=0
function dy=object1(t,y)
dy=zeros(4,1);
k1=2;k0=0.1;k2=0.1;k3=0.003; %k0 represents k-1
K=2; %K represents Michaelis constant Km,PDE
K1=0.0695;K4=0.4; %K1 represents KPDE
V1=1.26;m=0; %V1 represents VsGC,max，m represents orders
k4=K4*y(4)^m;
V2=K1*y(4);
dy(1)=-k1*y(1)*0.22+k0*y(3)+k4*y(2);
dy(3)=k1*y(1)*0.22-k0*y(3)-k2*y(3)-k3*y(3)*0.22;
dy(2)=k3*y(3)*0.22+k2*y(3)-k4*y(2);
dy(4)=V1*y(2)-y(4)*V2/(K+y(4)); %V2 represents VPDE,max
end

>> [T,Y]=ode45('object1',[0 30],[1 0 0 0]);
>> plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') %make images
>> xlabel('time(second)');ylabel('sGC forms');
>> [T,Y]=ode45('object1',[0 100],[1 0 0 0]);
>> gtext('D1');gtext('D2');gtext('D3');gtext('m=0');

>> V1=1.26;K1=0.0695;
>> Vp=V1*Y(:,2); %make the image of Vp-t
>> plot(Vp)
>> hold on
>> Vd=Y(:,4).^2*K1./(2+Y(:,4)); %make the image of Vd-t
>> plot(Vd)
>> xlabel('time(second)');ylabel('Vp and Vd(uM/s)');
>> gtext('Vp');gtext('Vd');gtext('m=0');

(2)When m=1
Just let m=1, V1=1.09, K1=0.032, K4=0.098 in (1).

(3)When m=2
Just let m=2, V1=0.8520, K1=0.0195, K4=0.011 in (1).

2.

function dy=object2(t,y)
dy=zeros(5,1);
k1=2;k0=0.1;k2=0.1;k3=0.003; %k0 represents k-1
K=2;L=0.22;K2=0.01; %K represents Michaelis constant Km,PDE
%L represents the initial quantity of NO in and outside cells
%K2 represents NO total consumption rate constant
K1=0.032;K4=0.098; %K1 represents KPDE
V1=1.09;m=1; %V1 represents VsGC,max，m represents orders
k4=K4*y(4)^m;
V2=K1*y(4);
dy(1)=-k1*y(1)*y(5)+k0*y(3)+k4*y(2);
dy(3)=k1*y(1)*y(5)-k0*y(3)-k2*y(3)-k3*y(3)*y(5);
dy(2)=k3*y(3)*y(5)+k2*y(3)-k4*y(2);
dy(4)=V1*y(2)-y(4)*V2/(K+y(4)); %V2 represents VPDE,max
dy(5)=L-K2*y(5);
end

>> [T,Y]=ode45('object2',[0 1000],[1 0 0 0 0.22]);
>> plot(Y(:,5),Y(:,4),'--') %make images of cGMP-NO
>> hold on
>> [T,Y]=ode45('object2',[0 1000],[1 0 0 0 0.22]);(***)
>> plot(Y(:,5),Y(:,4),'-')
>> xlabel('NO(uM)');ylabel('cGMP(uM)');
>> legend('no feedback','with feedback')

>> [T,Y]=ode45('object2',[0 1000],[1 0 0 0 0.22]);
>> plot(Y(:,5),Y(:,2),'--') %make images of sGC-NO
>> hold on
>> [T,Y]=ode45('object2',[0 1000],[1 0 0 0 0.22]);(***)
>> plot(Y(:,5),Y(:,2),'-')
>> xlabel('NO(uM)');ylabel('sGC(uM)');
>> legend('no feedback','with feedback')

Note

(***)represents when using this statement,we should change some parameters in the value in the M file, that is, let m=1, V1=1.09, K1=0.032, K4=0.098.

3.

function y=curve(s,t)
y=s(1)+s(2)*exp(s(3)*t);
end

% fit data
>> x=[185.5961,189.632,221.0417,287.9465,306.8128]; %experimental data
>> y=[11.78422893,12.51607142,13.93817204,17.98381967,20.11388292];
>> plot(x,y)
>> xlabel('sGC(U/L)');ylabel('cGMP(nmol/L)');
>> title('cGMP-sGC');

>> x=[185.5961,189.632,221.0417,287.9465,306.8128]; %experimental data
>> y=[11.78422893,12.51607142,13.93817204,17.98381967,20.11388292];
>> s0=[0.2 0.05 0.05]; %iteration initial value
>> sfit=lsqcurvefit('curve',s0,x,y); %least square curve fit
>> f=curve(sfit,x);
>> disp(sfit);

>> x=[185.5961,189.632,221.0417,287.9465,306.8128]; %experimental data
>> y=[11.78422893,12.51607142,13.93817204,17.98381967,20.11388292];
>> t=0:1:350;
>> a=8.3573;b=0.6748;k=0.0093;
>> yy=a+b*exp(k*t);
>> plot(x,y,'rx',t,yy,'g');
>> xlabel('sGC(U/L)');
>> ylabel('cGMP(nmol/L)');
>> title('cGMP-sGC');

4.

function dy=object3(t,y,Z)
lambda=10;s=0.000024;k=1;m=0.002;a=0.26;d=0.01;q=15;r=0.1;u=0.1;
dy=[lambda-d*y(1)-s*y(1)*y(3)/(1+r*y(3));
s*exp(-m*q)*Z(1,1)*Z(3,3)/(1+r*Z(3,3))-a*y(2);
k-u*y(3);];
end

>> lags=[15,1e-9,15];
>> history=[1;0;0];
>> tspan=[0,600];
>> sol=dde23(@object3,lags,history,tspan);
>> plot(sol.x,sol.y)

Note:repeated change the value of q, let q=15,5,2.5,1,0.1, and repeat the above code.

5.

function dy=object4(t,y)
dy=zeros(11,1);
V1=4;V2=2;c1=2;c2=2;
w1=2;w6=2;w7=2;w8=2;w9=2;w10=2;
k6=0.15;k7=0.15;k8=0.15;k9=0.15;k10=0.15;
m1=1;m2=2.5;m3=0.0103*1000;m4=0.0212*1000;m11=0.0171*1000;
d1=0.25;d2=0.25;d3=1.78;d4=1.78;d5=1.78;d11=0.1;
b6=0.5*7.5e-6;b7=0.4*7.5e-6;b8=0.3*7.5e-6;b9=0.2*7.5e-6;b10=0.1*7.5e-6;
dy(1)=(V1*y(5)/(m1+y(5)))*(y(11)/(c1+y(11)))-(b6*y(6)+b7*y(7)+b8*y(8)+b9*y(9)+b10*y(10))*y(1)-d1*y(1);
dy(2)=(V2*y(5))/(m2+y(5))*(y(11)/(c2+y(11)))-d2*y(2);
dy(3)=m11*y(11)-d3*y(3);
dy(4)=m3*y(3)-d4*y(4);
dy(5)=m4*y(4)-d5*y(5);
dy(6)=w6-k6*y(6);
dy(7)=w7-k7*y(7);
dy(8)=w8-k8*y(8);
dy(9)=w9-k9*y(9);
dy(10)=w10-k10*y(10);
dy(11)=w1-d11*y(11);
end

>> [T,Y]=ode45('object4',[0 100],[5 0 0 0 0 1 1 1 1 1 0]);
>> subplot(3,2,1);
>> plot(T,Y(:,1))
>> title('x1')
>> subplot(3,2,2);
>> plot(T,Y(:,2))
>> title('x2')
>> subplot(3,2,3);
>> plot(T,Y(:,3))
>> title('x3')
>> subplot(3,2,4);
>> plot(T,Y(:,4))
>> title('x4')
>> subplot(3,2,5);
>> plot(T,Y(:,5))
>> title('x5')

Team:SCUT-China/Modeling