Backup: Using Maxima to simulate a cell differentiation model 02

Two previous exercises are posted here as a record of my computational practice. 

 

Exercise with Maxima 3: simulating increase of k

 

In Huang's equation (Huang et al. 2007), assuming a1=a2=1, b1=b2=1, k1=k2=1, n=4 and all θs equal to 0.5, there will be three points where dx1/dt=dx2/dt=0, which are called "attractors" as stable states. They are [x1≈1.9961, x2≈0.0038], [x1=1.0, x2=1.0], and [x1≈0.0038, x2≈1.9961].

Here k is changed, increasing exponentially from 1.0 to 1.5, using the function k=1.5-0.5*e^(-0.005*t).

Here still ten replicates were conducted, with the starting states within the range [1.0±0.1, 1.0±0.1]. Same critieria as in last practice were used for determining the coefficients, except that the number of time intervals was set to be 1500 in order to let the trajectories reach the final attractors.

Figure. Trajectories generated for k increasing

The result shows that the trajectories here undergo a decrease in values of both x1 and x2 through the diagonal line at the first stage, before 'sliding' to the final attractors with lower x values (between 1.3 and 1.4) than those shown in last practice (higher than 1.8). Another thing is that some of the trajectories have tiny 'hooks' at the very beginning around the point [1.0, 1.0], which may also serve as a distinct feature.

Coding in Maxima:

f(u,v):= a*u^n/(th^n+u^n)+b*th^n/(th^n+v^n)-k*u;
a:1.0;
b:1.0;
k:1.0;
n:4;
th:0.5;
X1[0]:0.9+random(0.2);
Y1[0]:0.9+random(0.2);
X2[0]:0.9+random(0.2);
Y2[0]:0.9+random(0.2);
X3[0]:0.9+random(0.2);
Y3[0]:0.9+random(0.2);
X4[0]:0.9+random(0.2);
Y4[0]:0.9+random(0.2);
X5[0]:0.9+random(0.2);
Y5[0]:0.9+random(0.2);
X6[0]:0.9+random(0.2);
Y6[0]:0.9+random(0.2);
X7[0]:0.9+random(0.2);
Y7[0]:0.9+random(0.2);
X8[0]:0.9+random(0.2);
Y8[0]:0.9+random(0.2);
X9[0]:0.9+random(0.2);
Y9[0]:0.9+random(0.2);
X0[0]:0.9+random(0.2);
Y0[0]:0.9+random(0.2);
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y1[i]:Y1[i-1]+0.01*f(Y1[i-1],X1[i-1]),
      X1[i]:X1[i-1]+0.01*f(X1[i-1],Y1[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y2[i]:Y2[i-1]+0.01*f(Y2[i-1],X2[i-1]),
      X2[i]:X2[i-1]+0.01*f(X2[i-1],Y2[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y3[i]:Y3[i-1]+0.01*f(Y3[i-1],X3[i-1]),
      X3[i]:X3[i-1]+0.01*f(X3[i-1],Y3[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y4[i]:Y4[i-1]+0.01*f(Y4[i-1],X4[i-1]),
      X4[i]:X4[i-1]+0.01*f(X4[i-1],Y4[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y5[i]:Y5[i-1]+0.01*f(Y5[i-1],X5[i-1]),
      X5[i]:X5[i-1]+0.01*f(X5[i-1],Y5[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y6[i]:Y6[i-1]+0.01*f(Y6[i-1],X6[i-1]),
      X6[i]:X6[i-1]+0.01*f(X6[i-1],Y6[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y7[i]:Y7[i-1]+0.01*f(Y7[i-1],X7[i-1]),
      X7[i]:X7[i-1]+0.01*f(X7[i-1],Y7[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y8[i]:Y8[i-1]+0.01*f(Y8[i-1],X8[i-1]),
      X8[i]:X8[i-1]+0.01*f(X8[i-1],Y8[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y9[i]:Y9[i-1]+0.01*f(Y9[i-1],X9[i-1]),
      X9[i]:X9[i-1]+0.01*f(X9[i-1],Y9[i-1])];
k:1.0;
for i:1 while i<=1500
       do [
        k:1.5-0.5*%e^(-0.005*i),
      Y0[i]:Y0[i-1]+0.01*f(Y0[i-1],X0[i-1]),
      X0[i]:X0[i-1]+0.01*f(X0[i-1],Y0[i-1])];
x1:listarray(X1)$
y1:listarray(Y1)$
x2:listarray(X2)$
y2:listarray(Y2)$
x3:listarray(X3)$
y3:listarray(Y3)$
x4:listarray(X4)$
y4:listarray(Y4)$
x5:listarray(X5)$
y5:listarray(Y5)$
x6:listarray(X6)$
y6:listarray(Y6)$
x7:listarray(X7)$
y7:listarray(Y7)$
x8:listarray(X8)$
y8:listarray(Y8)$
x9:listarray(X9)$
y9:listarray(Y9)$
x0:listarray(X0)$
y0:listarray(Y0)$
plot2d([[discrete,x1,y1],[discrete,x2,y2],[discrete,x3,y3],[discrete,x4,y4],
[discrete,x5,y5],[discrete,x6,y6],[discrete,x7,y7],[discrete,x8,y8],
[discrete,x9,y9],[discrete,x0,y0]],[style,lines],[legend,"1","2","3","4",
"5","6","7","8","9","0"],[xlabel,"x1"],[ylabel,"x2"],
[gnuplot_term,png],[gnuplot_out_file,"graph1.png"]);