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]); |
![]() |
|
![]() |
|
![]() ![]() |
(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]); |
![]() ![]() |
|
![]() |
(2) |
Determine location of transcritical bifurcation by checking endemic equilibrium at J=0
> | endeq0:=simplify(eval(endeq,J=0)); |
![]() |
(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)); |
![]() |
(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']); |
![]() |
|
![]() |
> | 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']); |
![]() |
|
![]() |
|
![]() |
|
![]() |
> |
> | The Jacobian in the endemic equilbirium is
endJac:=simplify(eval(Jacobian(sys,[J,lII,lSS]),endeqsol)); |
![]() ![]() ![]() ![]() |
(5) |
> | Its characteristic polynomial is
endJacp:=CharacteristicPolynomial(endJac,lambda); |
![]() ![]() ![]() |
(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})); |
![]() |
|
![]() |
|
![]() |
|
![]() |
(7) |
> |