Stability analysis of the singular points and Hopf bifurcations of a tumor growth control model

Dániel András Drexler, Ilona Nagy, Valery G. Romanovski
Last modified: 01 October, 2023

1. a) Qualitative analysis

The tumor growth model

The following model describes the time evolution of the living tumor volume ((Subscript[x, 1])), the dead tumor volume ((Subscript[x, 2])), the drug level in the central compartment ((Subscript[x, 3])) and the drug level in the peripheral compartment ((Subscript[x, 4])):

(Subscript[x, 1])'=a (Subscript[x, 1])-n (Subscript[x, 1])-b (Subscript[x, 1]) ((Subscript[x, 3])/((Subscript[ED, 50])+(Subscript[x, 3])))
(Subscript[x, 2])'=n (Subscript[x, 1])-w (Subscript[x, 2])+b (Subscript[x, 1]) ((Subscript[x, 3])/((Subscript[ED, 50])+(Subscript[x, 3])))
(Subscript[x, 3])'=-c (Subscript[x, 3])-(Subscript[k, 1])(Subscript[x, 3])+(Subscript[k, 2])(Subscript[x, 4])+k((Subscript[x, 1])+(Subscript[x, 2]))
(Subscript[x, 4])'=(Subscript[k, 1])(Subscript[x, 3])-(Subscript[k, 2])(Subscript[x, 4])

The measured values in Table 2 are the following:

{a,b,n,w,ED,c,k1,k2}={(3131/12500),(3771/2500),(86643/1000000),(2993/100000),(17461/20000),(21471/25000),(8323/50),(8369/500)};
%//N

Calculations

System (1)

Quit

f1=a x1-n x1-b x1(x3/(ED+x3));
f2=n x1-w x2+b x1(x3/(ED+x3));
f3=-c x3-k1 x3+k2 x4+k(x1+x2);
f4=k1 x3- k2 x4;

Singular point (3)

sol=Solve[{f1,f2,f3,f4}==0,{x1,x2,x3,x4}]//Factor

The Jacobian at the nontrivial singular point (4), (5)

jac=D[{f1,f2,f3,f4},{{x1,x2,x3,x4}}]//FullSimplify;
{x1,x2,x3,x4}={x1,x2,x3,x4}/.sol[[2]];
jac=jac//FullSimplify
jac//MatrixForm

The characteristic polynomial (6)

pol=CharacteristicPolynomial[jac,y]//Factor

Variables[pol]

Initial parameter settings and characteristic polynomial (7)

{w,ED,c,k1,k2}={(2993/100000),(17461/20000),(21471/25000),(8323/50),(8369/500)};

pol=pol//Factor

The method to find pure imaginary eigenvalues

pb=(b0+b1 y+b2 (y^2))((y^2)+W)//Expand
eq=pol-pb;
uu=Series[eq,{y,0,7}]//FullSimplify
sys=LogicalExpand[uu==0]//FullSimplify

Elimination of (Subscript[b, 0]),(Subscript[b, 2]) in (9), (10)

b2=1;
eqs=Eliminate[sys,{b0,b1}]//FullSimplify

Setting the values a,b,n

ClearAll[a]; a=(3131/12500);
CylindricalDecomposition[eqs && x1>0 && x2>0 && x3>0 && x4>0 && n>0 && b>0 && W>0 &&a>n&&n>a-b, {b,n,k,W}]//FullSimplify

Reduce[((-140889196154970607369+7 (Sqrt[405097265038172811166168519191770258689]))/3213135150000)≤b<((9803161 (140889199374274976857+7 (Sqrt[405097265038172811166168519191770258689])))/6242123551227977106464450000),b];
%//N

Reduce[140889196154970607369+3213135150000 b>7 (Sqrt[405097265038172811166168519191770258689]),b];
%//N

The measured value is b=1.5084`.

ClearAll[b]; b=(3/2);
CylindricalDecomposition[eqs && x1>0 && x2>0 && x3>0 && x4>0 && n>0 && W>0 &&a>n&&n>a-b, {n,k,W}]//FullSimplify

Reduce[12488+25000 n==125 (NumericalApproximation[Root, (RowDefault["126.", FEExpressionsNumericalApproximationElider]), Root[2251099114179302908385820-140889200974673332369 ((#1)^2)+107104505 ((#1)^4)&,3], ]),n]
%//N

The measured value is n=0.086643`.

ClearAll[n];
n=((-12488+125 (NumericalApproximation[Root, (RowDefault["126.", FEExpressionsNumericalApproximationElider]), Root[2251099114179302908385820-140889200974673332369 ((#1)^2)+107104505 ((#1)^4)&,3], ]))/25000);

The eigenvalues (13)

Solve[pol==0,y]//N

1. b) Stability of the equilibrium state and the Hopf bifurcation

Step 1 - Shifting the singular point into the origin

The singular point of the system in the first quadrant is shifted into the origin and the parameters are chosen in such a way that the Jacobian at the origin has a pair of pure imaginary eigenvalues.

System (1)

Quit

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k,f1,f2,f3,f4];
f1=a x1-n x1-b x1(x3/(ED+x3));
f2=n x1-w x2+b x1(x3/(ED+x3));
f3=-c x3-k1 x3+k2 x4+k(x1+x2);
f4=k1 x3- k2 x4;

Singular points

sol=Solve[{f1==0,f2==0,f3==0,f4==0},{x1,x2,x3,x4}]//Factor
{y1,y2,y3,y4}={x1,x2,x3,x4}/.sol[[2]]

The singular point is shifted into the origin (17)

x1d=f1/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor
x2d=f2/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor
x3d=f3/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor
x4d=f4/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor

The Jacobian at the origin

Jac=D[{x1d,x2d,x3d,x4d},{{x1,x2,x3,x4}}]//Simplify;
JacOrigin=Jac/.{x1->0,x2->0,x3->0,x4->0}//Simplify

Eigenvalues of the Jacobian at the origin with the given parameter settings

{w,ED,c,k1,k2}={(2993/100000),(17461/20000),(21471/25000),(8323/50),(8369/500)};
a=(3131/12500);b=(3/2);n=((-12488+125 (NumericalApproximation[Root, (RowDefault["126.", FEExpressionsNumericalApproximationElider]), Root[2251099114179302908385820-140889200974673332369 ((#1)^2)+107104505 ((#1)^4)&,3], ]))/25000);
k=(1/100);
Eigenvalues[JacOrigin]//N

Step 2 - Taylor series expansion

Series expansion symbolically and the Jacobian at the origin

Quit

x1d=((a-b-n) (-a c ED w+c ED n w+(a^2) k x1-a b k x1-a k n x1+a k w x1-b k w x1-k n w x1) x3)/(k (a+w) (-b ED+a x3-b x3-n x3));
x2d=-(1/(k (a+w) (b ED-a x3+b x3+n x3)))(-(a^2) b ED k x1-a b ED k w x1+a b ED k w x2+b ED k (w^2) x2+(a^2) c ED w x3-a b c ED w x3-2 a c ED n w x3+b c ED n w x3+c ED (n^2) w x3+(a^2) b k x1 x3-a (b^2) k x1 x3+(a^2) k n x1 x3-2 a b k n x1 x3-a k (n^2) x1 x3+a b k w x1 x3-(b^2) k w x1 x3+a k n w x1 x3-2 b k n w x1 x3-k (n^2) w x1 x3-(a^2) k w x2 x3+a b k w x2 x3+a k n w x2 x3-a k (w^2) x2 x3+b k (w^2) x2 x3+k n (w^2) x2 x3);
x3d=k x1+k x2-c x3-k1 x3+k2 x4;
x4d=k1 x3-k2 x4;
r1=Series[x1d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal
r2=Series[x2d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal
r3=Series[x3d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal
r4=Series[x4d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal

Jac=D[{r1,r2,r3,r4},{{x1,x2,x3,x4}}]//Simplify;
JacOrigin=Jac/.{x1->0,x2->0,x3->0,x4->0}//Simplify

Series expansion with numerical values

A series expansion of the transformed system about the origin up to degree 4 will be used instead of the original system.

Quit

{w,ED,c,k1,k2}={(2993/100000),(17461/20000),(21471/25000),(8323/50),(8369/500)};
a=(3131/12500);b=(3/2);n=((-12488+125 (NumericalApproximation[Root, (RowDefault["126.", FEExpressionsNumericalApproximationElider]), Root[2251099114179302908385820-140889200974673332369 ((#1)^2)+107104505 ((#1)^4)&,3], ]))/25000);
k=(1/100);
x1d=((a-b-n) (-a c ED w+c ED n w+(a^2) k x1-a b k x1-a k n x1+a k w x1-b k w x1-k n w x1) x3)/(k (a+w) (-b ED+a x3-b x3-n x3));
x2d=-(1/(k (a+w) (b ED-a x3+b x3+n x3)))(-(a^2) b ED k x1-a b ED k w x1+a b ED k w x2+b ED k (w^2) x2+(a^2) c ED w x3-a b c ED w x3-2 a c ED n w x3+b c ED n w x3+c ED (n^2) w x3+(a^2) b k x1 x3-a (b^2) k x1 x3+(a^2) k n x1 x3-2 a b k n x1 x3-a k (n^2) x1 x3+a b k w x1 x3-(b^2) k w x1 x3+a k n w x1 x3-2 b k n w x1 x3-k (n^2) w x1 x3-(a^2) k w x2 x3+a b k w x2 x3+a k n w x2 x3-a k (w^2) x2 x3+b k (w^2) x2 x3+k n (w^2) x2 x3);
x3d=k x1+k x2-c x3-k1 x3+k2 x4;
x4d=k1 x3-k2 x4;
r1=Series[x1d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;
r2=Series[x2d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;
r3=Series[x3d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;
r4=Series[x4d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;

Separation of the linear and nonlinear parts

g1=r1/.{x1->t x1,x2->t x2,x3->t x3,x4->t x4};
rr1=Normal[Series[g1,{t,0,4}]]//Simplify;
g2=r2/.{x1->t x1,x2->t x2,x3->t x3,x4->t x4};
rr2=Normal[Series[g2,{t,0,4}]]//Simplify;
f1p=D[rr1,t]/.t->0//Simplify;
f2p=D[rr2,t]/.t->0//Simplify;
f1deg2=rr1-f1p t//Simplify//Factor;
f2deg2=rr2-f2p t//Simplify//Factor;

The full system

ff1=(SetPrecision[f1p,100] t+SetPrecision[f1deg2,100])/.t->1//Expand
ff2=(SetPrecision[f2p,100] t+SetPrecision[f2deg2,100])/.t->1//Expand
ff3=r3
ff4=r4

The full system with shorter numerical values in (18)

ff1//N
ff2//N
ff3//N
ff4//N

The nonlinear part

f1=SetPrecision[f1deg2,100]/.t->1//Expand
f2=SetPrecision[f2deg2,100]/.t->1//Expand
f3=0;
f4=0;

The Jacobian at the origin

Jac=D[{ff1,ff2,ff3,ff4},{{x1,x2,x3,x4}}]//Simplify;
JacOrigin=Jac/.{x1->0,x2->0,x3->0,x4->0}//Simplify;

SetPrecision[Eigenvalues[JacOrigin],100]

The Jacobian the singular point

J={{0,0,-((c (a-n) (-a+b+n) w)/(b k (a+w))),0},{a,-w,((c (a-n) (-a+b+n) w)/(b k (a+w))),0},{k,k,-c-k1,k2},{0,0,k1,-k2}}//FullSimplify;

SetPrecision[Eigenvalues[J],100]

Step 3 - Transformation to real Jordan normal form

The system is transformed in such a way that its linear part has real Jordan normal form (RJNF).
We remark that for the eigenvectors of J, both the exact values 
ev=Eigenvectors[J]//FullSimplify and the numerical values 
ev=SetPrecision[Eigenvectors[J],100] can be used. 
With SetPrecision the calculations are faster and provide more digits.

Transformation of J

ev=SetPrecision[Eigenvectors[J],100];

SetPrecision[Transpose[ev],100]//N//MatrixForm

The matrix S for which the product (S^(-1))J S has real Jordan normal form (19)

S=Transpose[{ev[[1]],ev[[2]],Re[ev[[3]]],Im[ev[[3]]]}];

SetPrecision[S,100]//N//MatrixForm

The inverse of S

Sinv=Inverse[S];

Verification that S transforms J to the real Jordan normal form (20)

RJNF=Sinv.J.S

SetPrecision[RJNF,100]

The entries (Subscript[a, 1]),(Subscript[a, 2]),β in (20)

a1=RJNF[[1,1]];
a2=RJNF[[2,2]];
beta=RJNF[[3,4]];

Transformation of the system

Consider the system (OverDot[x])=J x+f(x) where f contains the nonlinear terms.
With the substitution x=S y, this system is transformed into the form
(OverDot[y])=(S^(-1))J S y+(S^(-1))f(S y).
The product S y:

sy=S.{y1,y2,y3,y4}//FullSimplify;

The substitution f(S y):

FF1=f1/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify//Expand;
FF2=f2/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify//Expand;
FF3=f3/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify
FF4=f4/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify

The product (S^(-1))f(S y), that is, the nonlinear part in (OverDot[y]) in the new coordinates:

{h1,h2,h3,h4}=SetPrecision[Sinv.{FF1,FF2,FF3,FF4},100]//Simplify;

The full system in (21)

f1y=a1 y1+h1//Simplify;
f2y=a2 y2+h2//Simplify;
f3y=beta y4+h3//Simplify;
f4y=-beta y3+h4//Simplify;
jac=D[{f1y,f2y,f3y,f4y},{{y1,y2,y3,y4}}]/.{y1->0,y2->0,y3->0,y4->0};

jac//N//MatrixForm

Step 4 - Search for a Lyapunov function

We look for a polynomial   ϕ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]))=(Subscript[ϕ, 2])((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]))+(∑_(i+j+k+l=3)^4)(Subscript[p, ijkl])(Subscript[x, 1]^i)(Subscript[x, 2]^j)(Subscript[x, 3]^k)(Subscript[x, 4]^l)  
such that the quadratic part of ϕ,

(Subscript[ϕ, 2])((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]))=A (Subscript[x, 1]^2)+B (Subscript[x, 2]^2)+C (Subscript[x, 3]^2)+D (Subscript[x, 4]^2)

is positively defined and for the Lie derivative of the previous system we have that 

((∂ϕ)/(∂(Subscript[x, 1])))(Subscript[(OverDot[x]), 1])+((∂ϕ)/(∂(Subscript[x, 2])))(Subscript[(OverDot[x]), 2])+((∂ϕ)/(∂(Subscript[x, 3])))(Subscript[(OverDot[x]), 3])+((∂ϕ)/(∂(Subscript[x, 4])))(Subscript[(OverDot[x]), 4])=  (Subscript[g, 0]) [ ((Subscript[x, 1]^2)+(Subscript[x, 2]^2)) +    ((( (Subscript[x, 3]^2)+(Subscript[x, 4]^2)))^2)]+h.o.t.  

The following program compares the coefficients of the corresponding terms on both sides and thus the focus (or Lyapunov) quantity (Subscript[g, 0]) can be calculated.

The system in (21)

ClearAll[f1,f2,f3,f4];
f1=f1y/.{y1->x1,y2->x2,y3->x3,y4->x4};
f2=f2y/.{y1->x1,y2->x2,y3->x3,y4->x4};
f3=f3y/.{y1->x1,y2->x2,y3->x3,y4->x4};
f4=f4y/.{y1->x1,y2->x2,y3->x3,y4->x4};

Program

Ser[s_]:=Plus@@(Table[(x1^i)(x2^j)(x3^k) (x4^(s-i-j-k))p[i,j,k,s-i-j-k],{i,0,s},{j,0,s},{k,0,s},{l,0,s-i-j-k}]//Flatten//Union);
Hom[s_]:=Table[p[i,j,k,s-i-j-k],{i,0,s},{j,0,s},{k,0,s},{l,0,s-i-j-k}]//Flatten//Union;
hh=Sum[Ser[i],{i,3,4}];
V=hh+aa  (x1^2)+bb (x2^2)+ cc  (x3^2)+ dd  (x4^2);
Lie=D[V,x1] f1+D[V,x2] f2+D[V,x3] f3+D[V,x4] f4//Expand;
RHS=g0 (((x1^2)+(x2^2))+ ((((x3^2)+(x4^2)))^2))//Expand;
vv=Lie-RHS//Expand;
CoefPol[f_,s_]:=Module[{m,lis,t},lis={}; m=Expand[f];  Do[Do[Do[Do[If[i+j+k+l==s,lis=AppendTo[lis,Coefficient[m,(x1^i) (x2^j) (x3^k)(x4^l)]/.{x1->0,x2->0,x3->0,x4->0}]],{i,0,s}],{j,0,s}],{k,0,s}],{l,0,s}];
ls[s]=lis];
Do[CoefPol[vv,i],{i,1,4}]

Degree 2

ls[2]//FullSimplify

Variables[ls[2]]

sol2=Solve[ls[2]==0,{aa,bb,cc,dd}]//FullSimplify//Factor

{aa,bb,cc,dd}={aa,bb,cc,dd}/.sol2[[1]]

{aa,bb,cc,dd}//N

Positive definiteness of the quadratic part of V

S2=aa  (x1^2)+bb (x2^2)+ cc  (x3^2)+ cc (x4^2)

It can be seen that S2 will be positively defined if (Subscript[g, 0])<0 and c>0.

Degree 3

ls[3]//Factor;
sol3=Solve[ls[3]==0,Hom[3]]//Simplify;
{p[0,0,0,3],p[0,0,1,2],p[0,0,2,1],p[0,0,3,0],p[0,1,0,2],p[0,1,1,1],p[0,1,2,0],p[0,2,0,1],p[0,2,1,0],p[0,3,0,0],p[1,0,0,2],p[1,0,1,1],p[1,0,2,0],p[1,1,0,1],p[1,1,1,0],p[1,2,0,0],p[2,0,0,1],p[2,0,1,0],p[2,1,0,0],p[3,0,0,0]}={p[0,0,0,3],p[0,0,1,2],p[0,0,2,1],p[0,0,3,0],p[0,1,0,2],p[0,1,1,1],p[0,1,2,0],p[0,2,0,1],p[0,2,1,0],p[0,3,0,0],p[1,0,0,2],p[1,0,1,1],p[1,0,2,0],p[1,1,0,1],p[1,1,1,0],p[1,2,0,0],p[2,0,0,1],p[2,0,1,0],p[2,1,0,0],p[3,0,0,0]}/.sol3[[1]]//Simplify;
ls[3]//Simplify

Degree 4

ls[4]=ls[4]//Simplify;

sol4=Solve[ls[4]==0,AppendTo[Hom[4],g0]]

{g0,p[0,0,0,4],p[0,0,1,3],p[0,0,2,2],p[0,0,3,1],p[0,0,4,0],p[0,1,0,3],p[0,1,1,2],p[0,1,2,1],p[0,1,3,0],p[0,2,0,2],p[0,2,1,1],p[0,2,2,0],p[0,3,0,1],p[0,3,1,0],p[0,4,0,0],p[1,0,0,3],p[1,0,1,2],p[1,0,2,1],p[1,0,3,0],p[1,1,0,2],p[1,1,1,1],p[1,1,2,0],p[1,2,0,1],p[1,2,1,0],p[1,3,0,0],p[2,0,0,2],p[2,0,1,1],p[2,0,2,0],p[2,1,0,1],p[2,1,1,0],p[2,2,0,0],p[3,0,0,1],p[3,0,1,0],p[3,1,0,0],p[4,0,0,0]}={g0,p[0,0,0,4],p[0,0,1,3],p[0,0,2,2],p[0,0,3,1],p[0,0,4,0],p[0,1,0,3],p[0,1,1,2],p[0,1,2,1],p[0,1,3,0],p[0,2,0,2],p[0,2,1,1],p[0,2,2,0],p[0,3,0,1],p[0,3,1,0],p[0,4,0,0],p[1,0,0,3],p[1,0,1,2],p[1,0,2,1],p[1,0,3,0],p[1,1,0,2],p[1,1,1,1],p[1,1,2,0],p[1,2,0,1],p[1,2,1,0],p[1,3,0,0],p[2,0,0,2],p[2,0,1,1],p[2,0,2,0],p[2,1,0,1],p[2,1,1,0],p[2,2,0,0],p[3,0,0,1],p[3,0,1,0],p[3,1,0,0],p[4,0,0,0]}/.sol4[[1]];

g0

g0//N

From the above it can be seen that we can take e.g.  cc=1.

S2/.{ cc->1}

S2//N

V1=V  /.{ cc->1} //Simplify//Expand;

vd=D[V1,x1] f1+D[V1,x2] f2+D[V1,x3] f3+D[V1,x4] f4//Simplify;

gg=vd/.{x1->t x1,x2->t x2,x3->t x3,x4->t x4};
ff=Normal[Series[gg,{t,0,4}]]//Simplify;

MonomialList[ff,t]

S2/.{ cc->1}

From the above we can see that the lowest degree part of V is positively defined, that is, V is a positively defined Lyapunov function, whose Lie derivative is negatively defined.

1. c) Evolution of system states under perturbations

Plotting the original system

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots with the measured values

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k];
f1[x1_,x2_,x3_,x4_]:=a x1-n x1-b x1(x3/(ED+x3));
f2[x1_,x2_,x3_,x4_]:=n x1-w x2+b x1(x3/(ED+x3));
f3[x1_,x2_,x3_,x4_]:=-c x3-k1 x3+k2 x4+k(x1+x2);
f4[x1_,x2_,x3_,x4_]:=k1 x3- k2 x4;
{a,b,n,w,ED,c,k1,k2}={(3131/12500),(3771/2500),(86643/1000000),(2993/100000),(17461/20000),(21471/25000),(8323/50),(8369/500)};
k=(1/100);

ClearAll[nsol,ev,plotter3D];nsol=NSolve[Join@@Thread/@{{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]}==0,{x1,x2,x3,x4}>0},{x1,x2,x3,x4},20][[1]]
ev=Eigenvalues[D[{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]},{{x1,x2,x3,x4}}]/.nsol]
plotter3D[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,solution1,solution2,solution3,solution4,plot1,imagesize,imagesize2},
imagesize=350;imagesize2=250;
startingpoint=({x1,x2,x3,x4}/.nsol)+shift;
sys:=NDSolveValue[Join[{u1'[t]==f1[u1[t],u2[t],u3[t],u4[t]],u2'[t]==f2[u1[t],u2[t],u3[t],u4[t]],u3'[t]==f3[u1[t],u2[t],u3[t],u4[t]],u4'[t]==f4[u1[t],u2[t],u3[t],u4[t]]},Thread[{u1[0],u2[0],u3[0],u4[0]}==startingpoint]],
{u1,u2,u3,u4},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
solution4[t_]:=Delete[Through[sys[t]],4];
solution3[t_]:=Delete[Through[sys[t]],3];
solution2[t_]:=Delete[Through[sys[t]],2];
solution1[t_]:=Delete[Through[sys[t]],1];
{ParametricPlot3D[Evaluate[solution4[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution3[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution2[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution1[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 1]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 2]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[3]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[4]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2]}]

Trajectories going outwards, approaching the limit cycle from inside

plotter3D[1200,{0,0,0,0.5},Automatic,100,1000,1,Method->"BDF"]

Trajectories going inwards, approaching the limit cycle from outside

plotter3D[1000,{0,10,10,10},Automatic,100,1000,1,Method->"BDF"]

Positive perturbation

We perturb n by +1/50, so n≈0.152497. Each eigenvalue has a negative real part, the singular point is locally asymptotically stable.

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k];
f1[x1_,x2_,x3_,x4_]:=a x1-n x1-b x1(x3/(ED+x3));
f2[x1_,x2_,x3_,x4_]:=n x1-w x2+b x1(x3/(ED+x3));
f3[x1_,x2_,x3_,x4_]:=-c x3-k1 x3+k2 x4+k(x1+x2);
f4[x1_,x2_,x3_,x4_]:=k1 x3- k2 x4;
{w,ED,c,k1,k2}={(2993/100000),(17461/20000),(21471/25000),(8323/50),(8369/500)};
a=(3131/12500);b=(3/2);
n=((-12488+125 (NumericalApproximation[Root, (RowDefault["126.", FEExpressionsNumericalApproximationElider]), Root[2251099114179302908385820-140889200974673332369 ((#1)^2)+107104505 ((#1)^4)&,3], ]))/25000)+(1/50);
k=(1/100);

ClearAll[nsol,ev,plotter3D];nsol=NSolve[Join@@Thread/@{{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]}==0,{x1,x2,x3,x4}>0},{x1,x2,x3,x4},20][[1]]
ev=Eigenvalues[D[{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]},{{x1,x2,x3,x4}}]/.nsol]
plotter3D[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,solution1,solution2,solution3,solution4,plot1,imagesize,imagesize2},
imagesize=350;imagesize2=250;
startingpoint=({x1,x2,x3,x4}/.nsol)+shift;
sys:=NDSolveValue[Join[{u1'[t]==f1[u1[t],u2[t],u3[t],u4[t]],u2'[t]==f2[u1[t],u2[t],u3[t],u4[t]],u3'[t]==f3[u1[t],u2[t],u3[t],u4[t]],u4'[t]==f4[u1[t],u2[t],u3[t],u4[t]]},Thread[{u1[0],u2[0],u3[0],u4[0]}==startingpoint]],
{u1,u2,u3,u4},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
solution4[t_]:=Delete[Through[sys[t]],4];
solution3[t_]:=Delete[Through[sys[t]],3];
solution2[t_]:=Delete[Through[sys[t]],2];
solution1[t_]:=Delete[Through[sys[t]],1];
{ParametricPlot3D[Evaluate[solution4[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution3[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution2[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution1[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 1]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 2]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[3]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[4]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2]}]

Trajectories going inwards, towards the singular point

plotter3D[1000,{0,10,10,10},Automatic,100,1000,1,Method->"BDF"]

plotter3D[1200,{0,0,0,0.5},Automatic,100,1000,1,Method->"BDF"]

Negative perturbation

We perturb n by -1/50, so n≈0.112497. Two imaginary eigenvalues have positive real part, so the singular point becomes locally unstable and a stable limit cycle appears around it.

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k];
f1[x1_,x2_,x3_,x4_]:=a x1-n x1-b x1(x3/(ED+x3));
f2[x1_,x2_,x3_,x4_]:=n x1-w x2+b x1(x3/(ED+x3));
f3[x1_,x2_,x3_,x4_]:=-c x3-k1 x3+k2 x4+k(x1+x2);
f4[x1_,x2_,x3_,x4_]:=k1 x3- k2 x4;
{w,ED,c,k1,k2}={(2993/100000),(17461/20000),(21471/25000),(8323/50),(8369/500)};
a=(3131/12500);b=(3/2);
n=((-12488+125 (NumericalApproximation[Root, (RowDefault["126.", FEExpressionsNumericalApproximationElider]), Root[2251099114179302908385820-140889200974673332369 ((#1)^2)+107104505 ((#1)^4)&,3], ]))/25000)-(1/50);
k=(1/100);

ClearAll[nsol,ev,plotter3D];nsol=NSolve[Join@@Thread/@{{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]}==0,{x1,x2,x3,x4}>0},{x1,x2,x3,x4},20][[1]]
ev=Eigenvalues[D[{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]},{{x1,x2,x3,x4}}]/.nsol]
plotter3D[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,solution1,solution2,solution3,solution4,plot1,imagesize,imagesize2},
imagesize=350;imagesize2=250;
startingpoint=({x1,x2,x3,x4}/.nsol)+shift;
sys:=NDSolveValue[Join[{u1'[t]==f1[u1[t],u2[t],u3[t],u4[t]],u2'[t]==f2[u1[t],u2[t],u3[t],u4[t]],u3'[t]==f3[u1[t],u2[t],u3[t],u4[t]],u4'[t]==f4[u1[t],u2[t],u3[t],u4[t]]},Thread[{u1[0],u2[0],u3[0],u4[0]}==startingpoint]],
{u1,u2,u3,u4},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
solution4[t_]:=Delete[Through[sys[t]],4];
solution3[t_]:=Delete[Through[sys[t]],3];
solution2[t_]:=Delete[Through[sys[t]],2];
solution1[t_]:=Delete[Through[sys[t]],1];
{ParametricPlot3D[Evaluate[solution4[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution3[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution2[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution1[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 1]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 2]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[3]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[4]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2]}]

Case (i), Figures 1 and 3 (a)

The trajectories are going inward, approaching the limit cycle from outside.

plotter3D[1800,{0,10,10,10},Automatic,100,1000,1,Method->"BDF"]

Case (ii) , Figures 2 and 3 (b)

The trajectories are going outward, approaching the limit cycle from inside.

plotter3D[2500,{0,0,0,0.5},Automatic,100,1000,1,Method->"BDF"]

2. a) Qualitative analysis

The tumor growth model

Here we repeat the procedure described in the previous section.
The measured values in Table 3 are the following:

Quit

{a,b,n,w,ED,c,k1,k2}={(38437/100000),(22229/10000),(13881/100000),(71159/1000000),(117/80),(22711/100000),(21459/100),(13953/10)};
%//N

Calculations

System (1)

Quit

f1=a x1-n x1-b x1(x3/(ED+x3));
f2=n x1-w x2+b x1(x3/(ED+x3));
f3=-c x3-k1 x3+k2 x4+k(x1+x2);
f4=k1 x3- k2 x4;

Singular point (3)

sol=Solve[{f1,f2,f3,f4}==0,{x1,x2,x3,x4}]//Factor

The Jacobian at the nontrivial singular point (4), (5)

jac=D[{f1,f2,f3,f4},{{x1,x2,x3,x4}}]//FullSimplify;
{x1,x2,x3,x4}={x1,x2,x3,x4}/.sol[[2]];
jac=jac//FullSimplify
jac//MatrixForm

The characteristic polynomial (6)

pol=CharacteristicPolynomial[jac,y]//Factor

Variables[pol]

Initial parameter settings and characteristic polynomial

{w,ED,c,k1,k2}={(71159/1000000),(117/80),(22711/100000),(21459/100),(13953/10)};

pol=pol//Factor

The method to find pure imaginary eigenvalues

pb=(b0+b1 y+b2 (y^2))((y^2)+W)//Expand
eq=pol-pb;
uu=Series[eq,{y,0,7}]//FullSimplify
sys=LogicalExpand[uu==0]//FullSimplify

Elimination of (Subscript[b, 0]),(Subscript[b, 2])

b2=1;
eqs=Eliminate[sys,{b0,b1}]//FullSimplify

Setting the values a,b,n

ClearAll[a]; a=(38437/100000);
CylindricalDecomposition[eqs && x1>0 && x2>0 && x3>0 && x4>0 && n>0 && b>0 && W>0 && a>n && n>a-b, {b,n,k,W}]//FullSimplify

Reduce[17179639562122700701308600627400000 b>30151081 (361690835866668617225286079+1610188269 (Sqrt[50457031205553435413138389202156401])),b];
%//N

The measured value is b=2.2229, this is set for b.

ClearAll[b]; b=(22229/10000);
CylindricalDecomposition[eqs && x1>0 && x2>0 && x3>0 && x4>0 && n>0 && W>0 && a>n && n>a-b, {n,k,W}]//FullSimplify

Reduce[18177+25000 n==(NumericalApproximation[Root, (RowDefault["2.00", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "4"])]), Root[6427848213249358139332385841020380488175-16080051180984976488185730918002 ((#1)^2)+51714945568 ((#1)^4)&,3], ]),n]
%//N

The measured value is n=0.13881.

ClearAll[n];
n=((-18177+(NumericalApproximation[Root, (RowDefault["2.00", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "4"])]), Root[6427848213249358139332385841020380488175-16080051180984976488185730918002 ((#1)^2)+51714945568 ((#1)^4)&,3], ]))/25000);

The eigenvalues

Solve[pol==0,y]//N

2. b) Stability of the equilibrium state and the Hopf bifurcation

Step 1 - Shifting the singular point into the origin

The singular point of the system in the first quadrant is shifted into the origin and the parameters are chosen in such a way that the Jacobian at the origin has a pair of pure imaginary eigenvalues.

System (1)

Quit

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k,f1,f2,f3,f4];
f1=a x1-n x1-b x1(x3/(ED+x3));
f2=n x1-w x2+b x1(x3/(ED+x3));
f3=-c x3-k1 x3+k2 x4+k(x1+x2);
f4=k1 x3- k2 x4;

Singular points

sol=Solve[{f1==0,f2==0,f3==0,f4==0},{x1,x2,x3,x4}]//Factor
{y1,y2,y3,y4}={x1,x2,x3,x4}/.sol[[2]]

The singular point is shifted into the origin (17)

x1d=f1/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor
x2d=f2/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor
x3d=f3/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor
x4d=f4/.{x1->x1+y1,x2->x2+y2,x3->x3+y3,x4->x4+y4}//Factor

The Jacobian at the origin

Jac=D[{x1d,x2d,x3d,x4d},{{x1,x2,x3,x4}}]//Simplify;
JacOrigin=Jac/.{x1->0,x2->0,x3->0,x4->0}//Simplify

Eigenvalues of the Jacobian at the origin with the given parameter settings

{w,ED,c,k1,k2}={(71159/1000000),(117/80),(22711/100000),(21459/100),(13953/10)};
a=(38437/100000);b=(22229/10000);n=((-18177+(NumericalApproximation[Root, (RowDefault["2.00", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "4"])]), Root[6427848213249358139332385841020380488175-16080051180984976488185730918002 ((#1)^2)+51714945568 ((#1)^4)&,3], ]))/25000);
k=(1/100);
Eigenvalues[JacOrigin]//N

Step 2 - Taylor series expansion

Series expansion symbolically and the Jacobian at the origin

Quit

x1d=((a-b-n) (-a c ED w+c ED n w+(a^2) k x1-a b k x1-a k n x1+a k w x1-b k w x1-k n w x1) x3)/(k (a+w) (-b ED+a x3-b x3-n x3));
x2d=-(1/(k (a+w) (b ED-a x3+b x3+n x3)))(-(a^2) b ED k x1-a b ED k w x1+a b ED k w x2+b ED k (w^2) x2+(a^2) c ED w x3-a b c ED w x3-2 a c ED n w x3+b c ED n w x3+c ED (n^2) w x3+(a^2) b k x1 x3-a (b^2) k x1 x3+(a^2) k n x1 x3-2 a b k n x1 x3-a k (n^2) x1 x3+a b k w x1 x3-(b^2) k w x1 x3+a k n w x1 x3-2 b k n w x1 x3-k (n^2) w x1 x3-(a^2) k w x2 x3+a b k w x2 x3+a k n w x2 x3-a k (w^2) x2 x3+b k (w^2) x2 x3+k n (w^2) x2 x3);
x3d=k x1+k x2-c x3-k1 x3+k2 x4;
x4d=k1 x3-k2 x4;
r1=Series[x1d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal
r2=Series[x2d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal
r3=Series[x3d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal
r4=Series[x4d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal

Jac=D[{r1,r2,r3,r4},{{x1,x2,x3,x4}}]//Simplify;
JacOrigin=Jac/.{x1->0,x2->0,x3->0,x4->0}//Simplify

Series expansion with numerical values

A series expansion of the transformed system about the origin up to degree 4 will be used instead of the original system.

Quit

{w,ED,c,k1,k2}={(71159/1000000),(117/80),(22711/100000),(21459/100),(13953/10)};
a=(38437/100000);b=(22229/10000);n=((-18177+(NumericalApproximation[Root, (RowDefault["2.00", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "4"])]), Root[6427848213249358139332385841020380488175-16080051180984976488185730918002 ((#1)^2)+51714945568 ((#1)^4)&,3], ]))/25000);
k=(1/100);
x1d=((a-b-n) (-a c ED w+c ED n w+(a^2) k x1-a b k x1-a k n x1+a k w x1-b k w x1-k n w x1) x3)/(k (a+w) (-b ED+a x3-b x3-n x3));
x2d=-(1/(k (a+w) (b ED-a x3+b x3+n x3)))(-(a^2) b ED k x1-a b ED k w x1+a b ED k w x2+b ED k (w^2) x2+(a^2) c ED w x3-a b c ED w x3-2 a c ED n w x3+b c ED n w x3+c ED (n^2) w x3+(a^2) b k x1 x3-a (b^2) k x1 x3+(a^2) k n x1 x3-2 a b k n x1 x3-a k (n^2) x1 x3+a b k w x1 x3-(b^2) k w x1 x3+a k n w x1 x3-2 b k n w x1 x3-k (n^2) w x1 x3-(a^2) k w x2 x3+a b k w x2 x3+a k n w x2 x3-a k (w^2) x2 x3+b k (w^2) x2 x3+k n (w^2) x2 x3);
x3d=k x1+k x2-c x3-k1 x3+k2 x4;
x4d=k1 x3-k2 x4;
r1=Series[x1d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;
r2=Series[x2d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;
r3=Series[x3d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;
r4=Series[x4d,{x1,0,4},{x2,0,4},{x3,0,4},{x4,0,4}]//Normal;

Separation of the linear and nonlinear parts

g1=r1/.{x1->t x1,x2->t x2,x3->t x3,x4->t x4};
rr1=Normal[Series[g1,{t,0,4}]]//Simplify;
g2=r2/.{x1->t x1,x2->t x2,x3->t x3,x4->t x4};
rr2=Normal[Series[g2,{t,0,4}]]//Simplify;
f1p=D[rr1,t]/.t->0//Simplify;
f2p=D[rr2,t]/.t->0//Simplify;
f1deg2=rr1-f1p t//Simplify//Factor;
f2deg2=rr2-f2p t//Simplify//Factor;

The full system

ff1=(SetPrecision[f1p,100] t+SetPrecision[f1deg2,100])/.t->1//Expand
ff2=(SetPrecision[f2p,100] t+SetPrecision[f2deg2,100])/.t->1//Expand
ff3=r3
ff4=r4

The full system with shorter numerical values

ff1//N
ff2//N
ff3//N
ff4//N

The nonlinear part

f1=SetPrecision[f1deg2,100]/.t->1//Expand
f2=SetPrecision[f2deg2,100]/.t->1//Expand
f3=0;
f4=0;

The Jacobian at the origin

Jac=D[{ff1,ff2,ff3,ff4},{{x1,x2,x3,x4}}]//Simplify;
JacOrigin=Jac/.{x1->0,x2->0,x3->0,x4->0}//Simplify;

SetPrecision[Eigenvalues[JacOrigin],100]

The Jacobian the singular point

J={{0,0,-((c (a-n) (-a+b+n) w)/(b k (a+w))),0},{a,-w,((c (a-n) (-a+b+n) w)/(b k (a+w))),0},{k,k,-c-k1,k2},{0,0,k1,-k2}}//FullSimplify;

SetPrecision[Eigenvalues[J],100]

Step 3 - Transformation to real Jordan normal form

The system is transformed in such a way that its linear part has real Jordan normal form (RJNF).

Transformation of J

ev=SetPrecision[Eigenvectors[J],100];

SetPrecision[Transpose[ev],100]//N//MatrixForm

The matrix S for which the product (S^(-1))J S has real Jordan normal form

S=Transpose[{ev[[1]],ev[[2]],Re[ev[[3]]],Im[ev[[3]]]}];

SetPrecision[S,100]//N//MatrixForm

The inverse of S

Sinv=Inverse[S];

Verification that S transforms J to the real Jordan normal form

RJNF=Sinv.J.S

SetPrecision[RJNF,100]

The entries (Subscript[a, 1]),(Subscript[a, 2]),β

a1=RJNF[[1,1]];
a2=RJNF[[2,2]];
beta=RJNF[[3,4]];

Transformation of the system

Consider the system (OverDot[x])=J x+f(x) where f contains the nonlinear terms.
With the substitution x=S y, this system is transformed into the form
(OverDot[y])=(S^(-1))J S y+(S^(-1))f(S y).
The product S y:

sy=S.{y1,y2,y3,y4}//FullSimplify;

The substitution f(S y):

FF1=f1/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify//Expand;
FF2=f2/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify//Expand;
FF3=f3/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify
FF4=f4/.{x1->sy[[1]], x2->sy[[2]], x3->sy[[3]],x4->sy[[4]] }//Simplify

The product (S^(-1))f(S y), that is, the nonlinear part in (OverDot[y]) in the new coordinates:

{h1,h2,h3,h4}=SetPrecision[Sinv.{FF1,FF2,FF3,FF4},100]//Simplify;

The full system

f1y=a1 y1+h1//Simplify;
f2y=a2 y2+h2//Simplify;
f3y=beta y4+h3//Simplify;
f4y=-beta y3+h4//Simplify;
jac=D[{f1y,f2y,f3y,f4y},{{y1,y2,y3,y4}}]/.{y1->0,y2->0,y3->0,y4->0};

jac//N//MatrixForm

Step 4 - Search for a Lyapunov function

We look for a polynomial   ϕ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]))=(Subscript[ϕ, 2])((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]))+(∑_(i+j+k+l=3)^4)(Subscript[p, ijkl])(Subscript[x, 1]^i)(Subscript[x, 2]^j)(Subscript[x, 3]^k)(Subscript[x, 4]^l)  
such that the quadratic part of ϕ,

(Subscript[ϕ, 2])((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]))=A (Subscript[x, 1]^2)+B (Subscript[x, 2]^2)+C (Subscript[x, 3]^2)+D (Subscript[x, 4]^2)

is positively defined and for the Lie derivative of the previous system we have that 

((∂ϕ)/(∂(Subscript[x, 1])))(Subscript[(OverDot[x]), 1])+((∂ϕ)/(∂(Subscript[x, 2])))(Subscript[(OverDot[x]), 2])+((∂ϕ)/(∂(Subscript[x, 3])))(Subscript[(OverDot[x]), 3])+((∂ϕ)/(∂(Subscript[x, 4])))(Subscript[(OverDot[x]), 4])=  (Subscript[g, 0]) [ ((Subscript[x, 1]^2)+(Subscript[x, 2]^2)) +    ((( (Subscript[x, 3]^2)+(Subscript[x, 4]^2)))^2)]+h.o.t.  

The following program compares the coefficients of the corresponding terms on both sides and thus the focus (or Lyapunov) quantity (Subscript[g, 0]) can be calculated.

The system

ClearAll[f1,f2,f3,f4];
f1=f1y/.{y1->x1,y2->x2,y3->x3,y4->x4};
f2=f2y/.{y1->x1,y2->x2,y3->x3,y4->x4};
f3=f3y/.{y1->x1,y2->x2,y3->x3,y4->x4};
f4=f4y/.{y1->x1,y2->x2,y3->x3,y4->x4};

Program

Ser[s_]:=Plus@@(Table[(x1^i)(x2^j)(x3^k) (x4^(s-i-j-k))p[i,j,k,s-i-j-k],{i,0,s},{j,0,s},{k,0,s},{l,0,s-i-j-k}]//Flatten//Union);
Hom[s_]:=Table[p[i,j,k,s-i-j-k],{i,0,s},{j,0,s},{k,0,s},{l,0,s-i-j-k}]//Flatten//Union;
hh=Sum[Ser[i],{i,3,4}];
V=hh+aa  (x1^2)+bb (x2^2)+ cc  (x3^2)+ dd  (x4^2);
Lie=D[V,x1] f1+D[V,x2] f2+D[V,x3] f3+D[V,x4] f4//Expand;
RHS=g0 (((x1^2)+(x2^2))+ ((((x3^2)+(x4^2)))^2))//Expand;
vv=Lie-RHS//Expand;
CoefPol[f_,s_]:=Module[{m,lis,t},lis={}; m=Expand[f];  Do[Do[Do[Do[If[i+j+k+l==s,lis=AppendTo[lis,Coefficient[m,(x1^i) (x2^j) (x3^k)(x4^l)]/.{x1->0,x2->0,x3->0,x4->0}]],{i,0,s}],{j,0,s}],{k,0,s}],{l,0,s}];
ls[s]=lis];
Do[CoefPol[vv,i],{i,1,4}]

Degree 2

ls[2]//FullSimplify

Variables[ls[2]]

sol2=Solve[ls[2]==0,{aa,bb,cc,dd}]//FullSimplify//Factor

{aa,bb,cc,dd}={aa,bb,cc,dd}/.sol2[[1]]

{aa,bb,cc,dd}//N

Positive definiteness of the quadratic part of V

S2=aa  (x1^2)+bb (x2^2)+ cc  (x3^2)+ cc (x4^2)

It can be seen that S2 will be positively defined if (Subscript[g, 0])<0 and c>0.

Degree 3

ls[3]//Factor;
sol3=Solve[ls[3]==0,Hom[3]]//Simplify;
{p[0,0,0,3],p[0,0,1,2],p[0,0,2,1],p[0,0,3,0],p[0,1,0,2],p[0,1,1,1],p[0,1,2,0],p[0,2,0,1],p[0,2,1,0],p[0,3,0,0],p[1,0,0,2],p[1,0,1,1],p[1,0,2,0],p[1,1,0,1],p[1,1,1,0],p[1,2,0,0],p[2,0,0,1],p[2,0,1,0],p[2,1,0,0],p[3,0,0,0]}={p[0,0,0,3],p[0,0,1,2],p[0,0,2,1],p[0,0,3,0],p[0,1,0,2],p[0,1,1,1],p[0,1,2,0],p[0,2,0,1],p[0,2,1,0],p[0,3,0,0],p[1,0,0,2],p[1,0,1,1],p[1,0,2,0],p[1,1,0,1],p[1,1,1,0],p[1,2,0,0],p[2,0,0,1],p[2,0,1,0],p[2,1,0,0],p[3,0,0,0]}/.sol3[[1]]//Simplify;
ls[3]//Simplify

Degree 4

ls[4]=ls[4]//Simplify;

sol4=Solve[ls[4]==0,AppendTo[Hom[4],g0]]

{g0,p[0,0,0,4],p[0,0,1,3],p[0,0,2,2],p[0,0,3,1],p[0,0,4,0],p[0,1,0,3],p[0,1,1,2],p[0,1,2,1],p[0,1,3,0],p[0,2,0,2],p[0,2,1,1],p[0,2,2,0],p[0,3,0,1],p[0,3,1,0],p[0,4,0,0],p[1,0,0,3],p[1,0,1,2],p[1,0,2,1],p[1,0,3,0],p[1,1,0,2],p[1,1,1,1],p[1,1,2,0],p[1,2,0,1],p[1,2,1,0],p[1,3,0,0],p[2,0,0,2],p[2,0,1,1],p[2,0,2,0],p[2,1,0,1],p[2,1,1,0],p[2,2,0,0],p[3,0,0,1],p[3,0,1,0],p[3,1,0,0],p[4,0,0,0]}={g0,p[0,0,0,4],p[0,0,1,3],p[0,0,2,2],p[0,0,3,1],p[0,0,4,0],p[0,1,0,3],p[0,1,1,2],p[0,1,2,1],p[0,1,3,0],p[0,2,0,2],p[0,2,1,1],p[0,2,2,0],p[0,3,0,1],p[0,3,1,0],p[0,4,0,0],p[1,0,0,3],p[1,0,1,2],p[1,0,2,1],p[1,0,3,0],p[1,1,0,2],p[1,1,1,1],p[1,1,2,0],p[1,2,0,1],p[1,2,1,0],p[1,3,0,0],p[2,0,0,2],p[2,0,1,1],p[2,0,2,0],p[2,1,0,1],p[2,1,1,0],p[2,2,0,0],p[3,0,0,1],p[3,0,1,0],p[3,1,0,0],p[4,0,0,0]}/.sol4[[1]];

g0

g0//N

From the above it can be seen that we can take e.g.  cc=1.

S2/.{ cc->1}

S2//N

V1=V  /.{ cc->1} //Simplify//Expand;

vd=D[V1,x1] f1+D[V1,x2] f2+D[V1,x3] f3+D[V1,x4] f4//Simplify;

gg=vd/.{x1->t x1,x2->t x2,x3->t x3,x4->t x4};
ff=Normal[Series[gg,{t,0,4}]]//Simplify;

MonomialList[ff,t]

S2/.{ cc->1}

From the above we can see that the lowest degree part of V is positively defined, that is, V is a positively defined Lyapunov function, whose Lie derivative is negatively defined.

2. c) Evolution of system states under perturbations

Plotting the original system

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots with the measured values

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k];
f1[x1_,x2_,x3_,x4_]:=a x1-n x1-b x1(x3/(ED+x3));
f2[x1_,x2_,x3_,x4_]:=n x1-w x2+b x1(x3/(ED+x3));
f3[x1_,x2_,x3_,x4_]:=-c x3-k1 x3+k2 x4+k(x1+x2);
f4[x1_,x2_,x3_,x4_]:=k1 x3- k2 x4;
{a,b,n,w,ED,c,k1,k2}={(38437/100000),(22229/10000),(13881/100000),(71159/1000000),(117/80),(22711/100000),(21459/100),(13953/10)};
k=(1/100);
n//N

ClearAll[nsol,ev,plotter3D];nsol=NSolve[Join@@Thread/@{{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]}==0,{x1,x2,x3,x4}>0},{x1,x2,x3,x4},20][[1]]
ev=Eigenvalues[D[{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]},{{x1,x2,x3,x4}}]/.nsol]
plotter3D[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:100000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,solution1,solution2,solution3,solution4,plot1,imagesize,imagesize2},
imagesize=350;imagesize2=250;
startingpoint=({x1,x2,x3,x4}/.nsol)+shift;
sys:=NDSolveValue[Join[{u1'[t]==f1[u1[t],u2[t],u3[t],u4[t]],u2'[t]==f2[u1[t],u2[t],u3[t],u4[t]],u3'[t]==f3[u1[t],u2[t],u3[t],u4[t]],u4'[t]==f4[u1[t],u2[t],u3[t],u4[t]]},Thread[{u1[0],u2[0],u3[0],u4[0]}==startingpoint]],
{u1,u2,u3,u4},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
solution4[t_]:=Delete[Through[sys[t]],4];
solution3[t_]:=Delete[Through[sys[t]],3];
solution2[t_]:=Delete[Through[sys[t]],2];
solution1[t_]:=Delete[Through[sys[t]],1];
{ParametricPlot3D[Evaluate[solution4[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution3[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution2[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution1[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 1]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 2]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[3]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[4]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2]}]

Figure 4: Trajectories going inwards, converging to the singular point

plotter3D[1200,{0,0,1,1},Automatic,100,10000,1,Method->"BDF"]

Positive perturbation

If n is perturbed with a positive value then each eigenvalue has a negative real part  and the singular point is locally asymptotically stable.

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k];
f1[x1_,x2_,x3_,x4_]:=a x1-n x1-b x1(x3/(ED+x3));
f2[x1_,x2_,x3_,x4_]:=n x1-w x2+b x1(x3/(ED+x3));
f3[x1_,x2_,x3_,x4_]:=-c x3-k1 x3+k2 x4+k(x1+x2);
f4[x1_,x2_,x3_,x4_]:=k1 x3- k2 x4;
{w,ED,c,k1,k2}={(71159/1000000),(117/80),(22711/100000),(21459/100),(13953/10)};
a=(38437/100000);b=(22229/10000);
n=((-18177+(NumericalApproximation[Root, (RowDefault["2.00", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "4"])]), Root[6427848213249358139332385841020380488175-16080051180984976488185730918002 ((#1)^2)+51714945568 ((#1)^4)&,3], ]))/25000)+(1/100);
k=(1/100);
n//N

ClearAll[nsol,ev,plotter3D];nsol=NSolve[Join@@Thread/@{{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]}==0,{x1,x2,x3,x4}>0},{x1,x2,x3,x4},20][[1]]
ev=Eigenvalues[D[{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]},{{x1,x2,x3,x4}}]/.nsol]
plotter3D[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,solution1,solution2,solution3,solution4,plot1,imagesize,imagesize2},
imagesize=350;imagesize2=250;
startingpoint=({x1,x2,x3,x4}/.nsol)+shift;
sys:=NDSolveValue[Join[{u1'[t]==f1[u1[t],u2[t],u3[t],u4[t]],u2'[t]==f2[u1[t],u2[t],u3[t],u4[t]],u3'[t]==f3[u1[t],u2[t],u3[t],u4[t]],u4'[t]==f4[u1[t],u2[t],u3[t],u4[t]]},Thread[{u1[0],u2[0],u3[0],u4[0]}==startingpoint]],
{u1,u2,u3,u4},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
solution4[t_]:=Delete[Through[sys[t]],4];
solution3[t_]:=Delete[Through[sys[t]],3];
solution2[t_]:=Delete[Through[sys[t]],2];
solution1[t_]:=Delete[Through[sys[t]],1];
{ParametricPlot3D[Evaluate[solution4[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution3[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution2[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution1[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 1]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 2]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[3]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[4]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2]}]

Trajectories going inwards, towards the singular point

plotter3D[1000,{0,10,10,10},Automatic,100,1000,1,Method->"BDF"]

plotter3D[1200,{0,0,0.5,0.5},Automatic,100,1000,1,Method->"BDF"]

Negative perturbation

If n is perturbed with a negative value then two imaginary eigenvalues have positive real part, so the singular point becomes locally unstable and a stable limit cycle appears around it.

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[x1,x2,x3,a,b,n,w,ED,c,k1,k2,k];
f1[x1_,x2_,x3_,x4_]:=a x1-n x1-b x1(x3/(ED+x3));
f2[x1_,x2_,x3_,x4_]:=n x1-w x2+b x1(x3/(ED+x3));
f3[x1_,x2_,x3_,x4_]:=-c x3-k1 x3+k2 x4+k(x1+x2);
f4[x1_,x2_,x3_,x4_]:=k1 x3- k2 x4;
{w,ED,c,k1,k2}={(71159/1000000),(117/80),(22711/100000),(21459/100),(13953/10)};
a=(38437/100000);b=(22229/10000);
n=((-18177+(NumericalApproximation[Root, (RowDefault["2.00", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "4"])]), Root[6427848213249358139332385841020380488175-16080051180984976488185730918002 ((#1)^2)+51714945568 ((#1)^4)&,3], ]))/25000)-(5/100);
k=(1/100);
n//N

ClearAll[nsol,ev,plotter3D];nsol=NSolve[Join@@Thread/@{{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]}==0,{x1,x2,x3,x4}>0},{x1,x2,x3,x4},20][[1]]
ev=Eigenvalues[D[{f1[x1,x2,x3,x4],f2[x1,x2,x3,x4],f3[x1,x2,x3,x4],f4[x1,x2,x3,x4]},{{x1,x2,x3,x4}}]/.nsol]
plotter3D[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,solution1,solution2,solution3,solution4,plot1,imagesize,imagesize2},
imagesize=350;imagesize2=250;
startingpoint=({x1,x2,x3,x4}/.nsol)+shift;
sys:=NDSolveValue[Join[{u1'[t]==f1[u1[t],u2[t],u3[t],u4[t]],u2'[t]==f2[u1[t],u2[t],u3[t],u4[t]],u3'[t]==f3[u1[t],u2[t],u3[t],u4[t]],u4'[t]==f4[u1[t],u2[t],u3[t],u4[t]]},Thread[{u1[0],u2[0],u3[0],u4[0]}==startingpoint]],
{u1,u2,u3,u4},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
solution4[t_]:=Delete[Through[sys[t]],4];
solution3[t_]:=Delete[Through[sys[t]],3];
solution2[t_]:=Delete[Through[sys[t]],2];
solution1[t_]:=Delete[Through[sys[t]],1];
{ParametricPlot3D[Evaluate[solution4[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution3[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution2[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 1]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
ParametricPlot3D[Evaluate[solution1[t]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"(Subscript[x, 2]) [(mm^3)]","(Subscript[x, 3]) [mg/kg]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[14],ImageSize->imagesize,BoxRatios->{1,1,1}],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 1]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 2]) [(mm^3)]"},LabelStyle->Directive[12],ImagePadding->{{22,48},{10,22}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[3]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 3]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2],
Plot[Evaluate[solution[t][[4]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{"t [days]","(Subscript[x, 4]) [mg/kg]"},LabelStyle->Directive[12],ImagePadding->{{25,48},{10,20}},ImageSize->imagesize2]}]

Figure 5. (a)

The trajectories are going inward, approaching the limit cycle from outside.

plotter3D[1300,{0,1,1,1},Automatic,100,1000,1,Method->"BDF"]

Figure 5. (b)

The trajectories are going outward, approaching the limit cycle from inside.

plotter3D[1300,{0,0,0,0.2},Automatic,100,1000,1,Method->"BDF"]