SIS-ode.mw

ODE mean field model for SIS network 

> restart:with(LinearAlgebra):with(VectorCalculus):
 

RHS of ODE, using J fo rnumber of infected to distignuish from imaginary unit 

> S:=1-J;
lSI:=k/2-lSS-lII;
sys:=Vector([-r*J+p*lSI,p*lSI*(lSI/S+1)-2*r*lII,(r+w)*lSI-2*p*lSI*lSS/S]);
 

 

 

S := VectorCalculus:-`+`(1, VectorCalculus:-`-`(J))
lSI := VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(`/`(1, 2), k), VectorCalculus:-`-`(lSS)), VectorCalculus:-`-`(lII))
sys := VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`-`(VectorCalculus:-`*`(r, J)), VectorCalculus:-`*`(p, VectorCalculus:-`+`(VectorCalculus:-`+`(Ve...
sys := VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`-`(VectorCalculus:-`*`(r, J)), VectorCalculus:-`*`(p, VectorCalculus:-`+`(VectorCalculus:-`+`(Ve...
(1)
 

endeq are [LII,LSS,p] component of endemic equilibrium, parametrized by J 

> endeqsol:=solve({sys[1],sys[2],sys[3]},{p,lSS,lII}):
endeq:=subs(endeqsol,[lII,lSS,p]);factor(endeq[3]);
 

 

endeq := [VectorCalculus:-`-`(VectorCalculus:-`*`(`/`(1, 2), VectorCalculus:-`*`(VectorCalculus:-`*`(J, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-...
endeq := [VectorCalculus:-`-`(VectorCalculus:-`*`(`/`(1, 2), VectorCalculus:-`*`(VectorCalculus:-`*`(J, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-...
VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(r, VectorCalculus:-`*`(`*`(`^`(J, 2)), w)), w), VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(2, w), J))),... (2)
 

Determine location of transcritical bifurcation by checking endemic equilibrium at J=0 

> endeq0:=simplify(eval(endeq,J=0));
 

endeq0 := [0, VectorCalculus:-`*`(`/`(1, 2), k), VectorCalculus:-`*`(VectorCalculus:-`+`(r, w), `/`(1, `*`(k)))] (3)
 

> Check derivative of endemic equilibrium to check criticality of transcritical bifurcation (is endemic equilibrium stable or unstable?
simplify(eval(diff(endeq[3],J),J=0));
 

VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`-`(VectorCalculus:-`*`(k, w)), VectorCalculus:-`*`(k, r)), r), w), `/`(1, `*`(`^`(k, 2)))) (4)
 

> For the concrete values of k=10, w=0.04 and r=0.002 this is (giving the plot below)
pI:=eval(endeq[3],{r=0.002,w=0.04,k=10});
plot([pI,J,J=0..0.98],view=[0..0.008,0..1],labels=['p','I']);
 

 

pI := VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`+`(0.42e-1, VectorCalculus:-`*`(0.4e-1, `*`(`^`(J, 2)))), VectorCalculus:-`-`(VectorCalculus:-`*`(0.8e-1, J))), `/`(1, `*`(VectorCalculus...
Plot_2d
 

> The fold can be obtained by simply checking where the derivative of p wrt J is zero. We parametrize the fold curve in the (p,w) plane by J. The plot below is the fold curve in (p,w,J)-space.
wfold:=solve(diff(endeq[3],J)=0,w);
pfold:=simplify(eval(endeq[3],w=wfold));
fc:=eval([pfold,wfold],{k=10,r=0.002});
plots[spacecurve]([fc[1],fc[2],J,J=0..0.9],view=[0..0.003,0..0.05,0..1],axes=normal,labels=['p','w','I']);
 

 

 

 

wfold := VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(r, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`-`(k), -1), VectorCalculus:-`*`(2, J))), `/`(1, `*`(VectorCalculus:-`*`...
pfold := VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(2, r), `/`(1, `*`(VectorCalculus:-`*`(VectorCalculus:-`+`(-1, J), VectorCalculus:-`+`(-1, k))))))
fc := [VectorCalculus:-`-`(VectorCalculus:-`*`(0.4444444444e-3, `/`(1, `*`(VectorCalculus:-`+`(-1, J))))), VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(0.2222222222e-3, VectorCalculus:-...
Plot_2d
 

>
 

> The Jacobian in the endemic equilbirium is
endJac:=simplify(eval(Jacobian(sys,[J,lII,lSS]),endeqsol));
 

endJac := Matrix(%id = 5503597250)
endJac := Matrix(%id = 5503597250)
endJac := Matrix(%id = 5503597250)
endJac := Matrix(%id = 5503597250)
(5)
 

> Its characteristic polynomial is
endJacp:=CharacteristicPolynomial(endJac,lambda);

 

endJacp := VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(`*`(`^`(lambda, 3)), VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCa...
endJacp := VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(`*`(`^`(lambda, 3)), VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCa...
endJacp := VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(`*`(`^`(lambda, 3)), VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCa...
(6)
 

Gross et al (2006) observed a Hopf bifurcation above a certain minimal w. We only calculate this minimal w by checking for a double-eigenvalue zero of the Jacobian. At wdpzero dp/dlambda(0)=0. This has to be satisfied simultaneously with p(0)=0. The resulting equation has 3 solutions for J, only one of which is between 0 and 1. 

> wdpzero:=solve(eval(diff(endJacp,lambda),lambda=0),w);
solve(wfold=wdpzero,J);
J00:=[%][1];
wHopfmin:=evalf(eval(eval(wfold,J=J00),{k=10,r=0.002}));
 

 

 

 

wdpzero := VectorCalculus:-`*`(VectorCalculus:-`*`(r, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(J, VectorCalculus:-`*`(3, `*`(`^`(J, 2)))), VectorCalculus:-`-`(Ve...
VectorCalculus:-`*`(`/`(1, 3), VectorCalculus:-`*`(VectorCalculus:-`+`(-1, sqrt(VectorCalculus:-`+`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(VectorCalculus:-`*`(6, k))), VectorCalculus:-`*`(6, `*`(`...
J00 := VectorCalculus:-`*`(`/`(1, 3), VectorCalculus:-`*`(VectorCalculus:-`+`(-1, sqrt(VectorCalculus:-`+`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(VectorCalculus:-`*`(6, k))), VectorCalculus:-`*`(6...
wHopfmin := 0.6740853313e-1 (7)
 

>