Backup: Using Maxima to simulate a cell differentiation model 01

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

Exercise with Maxima 2: simulating increase of θ

 

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].

However it is clear that when θ increases from 0.5 to 1.0, the central point [1.0, 1.0] will not be an attractor anymore. If the initial state of a cell is near this point, it will slide away and finally fall into either of the two side attractors. So what will the trajectory be like?

Here, θ was set to increase from 0.5 to 1.0 at a constant rate over time. As the unit of time was unclarified in the original equation, proper value of the time was determined with the aim that the final trajectory can reach the final attractor without stopping on halfway. The coefficient between ∆x and ∆t is chosen cautiously to prevent x values from falling to the transient attractor nearby before t is changed again.

After each time interval, ∆t, θ increases by 0.0005 and the x1, x2 values change by ∆x=10*x'(t)*∆t so that discrete values of x1 and x2 were obtained to draw the trajectory. Ten replicates were conducted with the initial states randomly picked within the range [1.0±0.1, 1.0±0.1]. The simulation was performed with Maxima 5.24.0 (wxMaxima 11.04.0). The result is shown in the figure below and the source code is given at the bottom.

The conclusion is that when a cell's state goes from around the central state towards either one of the two side attractors due to the increase of θ, its trajectory goes towards to the central point [1.0, 1.0] shortly at the very beginning, before they are then attracted away. In a realistic experiment, the initial direction of the trajectory should just be at random because we don't know which way is towards the supposed central state.

Coding in Maxima:

f(u,v):= a*u^n/(th^n+u^n)+b*th^n/(th^n+v^n)-k*u;
a:1;
b:1;
k:1;
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<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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])];
th:0.5;
for i:1 while i<=1000
       do [
        th:th+0.0005,
      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"]);