SIS-models 

If you click on !!! (the menu bar at the top) maple will execute all commands.  

You can also enter each command as you go along  reading this document. You may enter your own commands at any time and place in this document. 

> restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1
 

> `:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
`:=`(iterate, proc (F, P0, t) option remember; description
 

> `:=`(timeplot2, proc (F, P0, T) local orb; description
`:=`(timeplot2, proc (F, P0, T) local orb; description
`:=`(timeplot2, proc (F, P0, T) local orb; description
`:=`(timeplot2, proc (F, P0, T) local orb; description
`:=`(timeplot2, proc (F, P0, T) local orb; description
`:=`(timeplot2, proc (F, P0, T) local orb; description
`:=`(timeplot2, proc (F, P0, T) local orb; description
`:=`(timeplot2, proc (F, P0, T) local orb; description
 

> `:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
 

> `:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
`:=`(dirfield, proc (F, Pmin, Pmax, Qmin, Qmax) description `(`+`(Pmin, `*`(0.5e-1, `*`(i))), `+`(Qmin, `*`(0.5e-1..." align="center" border="0">
 

> `:=`(phasedia2, proc (F, P0, T) local orb; description
`:=`(phasedia2, proc (F, P0, T) local orb; description
`:=`(phasedia2, proc (F, P0, T) local orb; description
`:=`(phasedia2, proc (F, P0, T) local orb; description
`:=`(phasedia2, proc (F, P0, T) local orb; description
`:=`(phasedia2, proc (F, P0, T) local orb; description
`:=`(phasedia2, proc (F, P0, T) local orb; description
`:=`(phasedia2, proc (F, P0, T) local orb; description
 

> `:=`(SIS, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`(alpha, `*`(P[1], `*`(P[2])))), `*`(gamma, `*`(P[2]))), `+`(P[2], `*`(alpha, `*`(..." align="center" border="0">
`:=`(SIS, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`(alpha, `*`(P[1], `*`(P[2])))), `*`(gamma, `*`(P[2]))), `+`(P[2], `*`(alpha, `*`(..." align="center" border="0">
`:=`(SIS, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`(alpha, `*`(P[1], `*`(P[2])))), `*`(gamma, `*`(P[2]))), `+`(P[2], `*`(alpha, `*`(..." align="center" border="0">
`:=`(SIS, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`(alpha, `*`(P[1], `*`(P[2])))), `*`(gamma, `*`(P[2]))), `+`(P[2], `*`(alpha, `*`(..." align="center" border="0">
`:=`(SIS, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`(alpha, `*`(P[1], `*`(P[2])))), `*`(gamma, `*`(P[2]))), `+`(P[2], `*`(alpha, `*`(..." align="center" border="0">
`:=`(SIS, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`(alpha, `*`(P[1], `*`(P[2])))), `*`(gamma, `*`(P[2]))), `+`(P[2], `*`(alpha, `*`(..." align="center" border="0">
 

> `:=`(param, proc (F) local a, b; description `(0, 1))[1]); `:=`(a, `+`(1, `-`(F(`<,>`(1, 1))[1]), b)); return [a, b] end proc); -1" align="center" border="0">
`:=`(param, proc (F) local a, b; description `(0, 1))[1]); `:=`(a, `+`(1, `-`(F(`<,>`(1, 1))[1]), b)); return [a, b] end proc); -1" align="center" border="0">
 

> param(SIS(alpha, gamma)); 1
 

[alpha, gamma] (1)
 

General theory 

> equi2(proc (P) options operator, arrow; `<,>`((SIS(alpha, gamma))(P)[1], (SIS(alpha, gamma))(P)[2]) end proc); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = P[1], P[2] = 0], [P[1] = `/`(`*`(gamma), `*`(alpha)), P[2] = P[2]]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 141068260), Vector[column](%id = 138873452)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 140954500), Vector[column](%id = 140547452) (2)
 

In a SIS-model, I = P[2] follows logistic growth because P[1]+P[2]=N is constant 

> subs(P[1] = VectorCalculus:-`+`(N, VectorCalculus:-`-`(P[2])), (SIS(alpha, gamma))(P)[2]); 1
 

`+`(P[2], `*`(alpha, `*`(`+`(N, `-`(P[2])), `*`(P[2]))), `-`(`*`(gamma, `*`(P[2])))) (3)
 

> % = VectorCalculus:-`+`(P[2], VectorCalculus:-`*`(VectorCalculus:-`*`(alpha, P[2]), VectorCalculus:-`+`(VectorCalculus:-`+`(N, VectorCalculus:-`-`(VectorCalculus:-`*`(gamma, `/`(1, `*`(alpha))))), Vec...
 

`+`(P[2], `*`(alpha, `*`(`+`(N, `-`(P[2])), `*`(P[2]))), `-`(`*`(gamma, `*`(P[2])))) = `+`(P[2], `*`(alpha, `*`(P[2], `*`(`+`(N, `-`(`/`(`*`(gamma), `*`(alpha))), `-`(P[2])))))) (4)
 

> VectorCalculus:-`+`(P[2], VectorCalculus:-`*`(VectorCalculus:-`*`(alpha, P[2]), VectorCalculus:-`+`(VectorCalculus:-`+`(N, VectorCalculus:-`-`(VectorCalculus:-`*`(gamma, `/`(1, `*`(alpha))))), VectorC...
 

`+`(P[2], `*`(alpha, `*`(P[2], `*`(`+`(N, `-`(`/`(`*`(gamma), `*`(alpha))), `-`(P[2])))))) = `+`(P[2], `*`(alpha, `*`(`+`(N, `-`(`/`(`*`(gamma), `*`(alpha)))), `*`(P[2](`+`(1, `-`(`/`(`*`(P[2]), `*`(`... (5)
 

This is logistic growth with carrying capacity N-γ/α or N-ρ where ρ=γ/α is the relative removal rate. 

Example 1 

> `:=`(F, SIS(0.1e-1, .25)); -1
 

Relative removal rate and contact number 

> `:=`(rho, VectorCalculus:-`*`(param(F)[2], `/`(1, `*`(param(F)[1])))); 1
 

25. (6)
 

> `:=`(sigma, VectorCalculus:-`*`(100, `/`(1, `*`(rho)))); 1
 

4.000000000 (7)
 

> display([timeplot2(F, `<,>`(95, 5), 25), plot(VectorCalculus:-`+`(100, VectorCalculus:-`-`(rho)), 0 .. 25, color = magenta)]); 1
 

Plot_2d
 

Example 2 

> `:=`(F, SIS(0.5e-2, .25)); -1
 

Relative removal rate and contact number 

> `:=`(rho, VectorCalculus:-`*`(param(F)[2], `/`(1, `*`(param(F)[1])))); 1
 

50.00000000 (8)
 

> `:=`(sigma, VectorCalculus:-`*`(100, `/`(1, `*`(rho)))); 1
 

2.000000000 (9)
 

> display([timeplot2(F, `<,>`(95, 5, 0), 25), plot(VectorCalculus:-`+`(100, VectorCalculus:-`-`(rho)), 0 .. 25, color = magenta)]); 1
 

Plot_2d
 

Example 3 

> `:=`(F, SIS(0.3e-2, .25)); -1
 

Relative removal rate and contact number 

> `:=`(rho, VectorCalculus:-`*`(param(F)[2], `/`(1, `*`(param(F)[1])))); 1
 

83.33333332 (10)
 

> `:=`(sigma, VectorCalculus:-`*`(100, `/`(1, `*`(rho)))); 1
 

1.200000000 (11)
 

> display([timeplot2(F, `<,>`(95, 5, 0), 75), plot(VectorCalculus:-`+`(100, VectorCalculus:-`-`(rho)), 0 .. 75, color = magenta)]); 1
 

Plot_2d
 

Differentiated SIS model 

> `:=`(SISto, proc (af, gf, Nf, am, gm, Nm) description `(`+`(P[1], `*`(af, `*`(`+`(Nf, `-`(P[1])), `*`(P[2]))), `-`(`*`(gf, `*`(..." align="center" border="0">
`:=`(SISto, proc (af, gf, Nf, am, gm, Nm) description `(`+`(P[1], `*`(af, `*`(`+`(Nf, `-`(P[1])), `*`(P[2]))), `-`(`*`(gf, `*`(..." align="center" border="0">
`:=`(SISto, proc (af, gf, Nf, am, gm, Nm) description `(`+`(P[1], `*`(af, `*`(`+`(Nf, `-`(P[1])), `*`(P[2]))), `-`(`*`(gf, `*`(..." align="center" border="0">
`:=`(SISto, proc (af, gf, Nf, am, gm, Nm) description `(`+`(P[1], `*`(af, `*`(`+`(Nf, `-`(P[1])), `*`(P[2]))), `-`(`*`(gf, `*`(..." align="center" border="0">
`:=`(SISto, proc (af, gf, Nf, am, gm, Nm) description `(`+`(P[1], `*`(af, `*`(`+`(Nf, `-`(P[1])), `*`(P[2]))), `-`(`*`(gf, `*`(..." align="center" border="0">
 

General theory 

> `:=`(F, SISto(alpha[f], gamma[f], N[f], alpha[m], gamma[m], N[m])); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(P[1], VectorCalculus:-`*`(VectorCalculus:-`*`(alpha[f], VectorCalculus:-`+`(N[f], VectorCalculus:-`-`(P[...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(P[1], VectorCalculus:-`*`(VectorCalculus:-`*`(alpha[f], VectorCalculus:-`+`(N[f], VectorCalculus:-`-`(P[...
(12)
 

> F(P); 1
 

Vector[column](%id = 140738728) (13)
 

The equilibria  

> simplify(solve({F(P)[1] = P[1], F(P)[2] = P[2]}, [P[1], P[2]])); 1
 

[[P[1] = 0, P[2] = 0], [P[1] = `/`(`*`(`+`(`*`(alpha[m], `*`(N[m], `*`(alpha[f], `*`(N[f])))), `-`(`*`(gamma[m], `*`(gamma[f]))))), `*`(alpha[m], `*`(`+`(`*`(N[m], `*`(alpha[f])), gamma[f])))), P[2] =... (14)
 

> equi2(proc (P) options operator, arrow; `<,>`(F(P)[1], F(P)[2]) end proc); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 0, P[2] = 0], [P[1] = `+`(`-`(`/`(`*`(`+`(`-`(`*`(alpha[m], `*`(N[m], `*`(alpha[f], `*`(N[f]))))), `*`(gamma[m], `*`(gamma[f])))), `*`(alpha[m], `*`(`+`(`*`(N[m], `*`(alpha[f])), gamma[f]))))...
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
Vector[column](%id = 138997116), Vector[column](%id = 135338164)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
Vector[column](%id = 137329256), Vector[column](%id = 141281052)
(15)
 

Example 1 

> `:=`(F, SISto(0.2e-2, 0.3e-1, 100, 0.4e-2, 0.1e-1, 100)); -1
 

The equilibria are 

> equi2(F); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 0., P[2] = 0.], [P[1] = 86.63043478, P[2] = 97.19512195]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 138835384), Vector[column](%id = 142345404)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 142339044), Vector[column](%id = 142339116) (16)
 

> timeplot2(F, `<,>`(30, 50), 20); 1
 

Plot_2d
 

> timeplot2(F, `<,>`(40, 1), 20); 1
 

Plot_2d
 

Notice that at the outset the disease is decreasing in one part of the population and increasing in the other part.   

> `:=`(df, dirfield(F, 86, 88, 96, 98)); -1
 

> `:=`(orb, seq(phasedia2(F, P0, 10), P0 = `union`(`union`(`union`({seq(`<,>`(VectorCalculus:-`+`(86, VectorCalculus:-`*`(.2, k)), 96), k = 0 .. 10)}, {seq(`<,>`(88, VectorCalculus:-`+`(96, VectorCalcul...
`:=`(orb, seq(phasedia2(F, P0, 10), P0 = `union`(`union`(`union`({seq(`<,>`(VectorCalculus:-`+`(86, VectorCalculus:-`*`(.2, k)), 96), k = 0 .. 10)}, {seq(`<,>`(88, VectorCalculus:-`+`(96, VectorCalcul...
 

> display(df); 1
 

Plot_2d
 

> display(orb); 1
 

Plot_2d
 

> display(df, orb); 1
 

Plot_2d
 

> display(seq(phasedia2(F, P0, 10), P0 = `union`(`union`(`union`({seq(`<,>`(VectorCalculus:-`+`(25, VectorCalculus:-`*`(2, k)), 15), k = 0 .. 45)}, {seq(`<,>`(25, VectorCalculus:-`+`(15, VectorCalculus:...
display(seq(phasedia2(F, P0, 10), P0 = `union`(`union`(`union`({seq(`<,>`(VectorCalculus:-`+`(25, VectorCalculus:-`*`(2, k)), 15), k = 0 .. 45)}, {seq(`<,>`(25, VectorCalculus:-`+`(15, VectorCalculus:...
 

Plot_2d
 

Example 2 

> `:=`(F, SISto(0.2e-1, 0.5e-1, 100, 0.2e-1, 0.5e-1, 100)); -1
 

The equlibria are 

> equi2(F); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 0., P[2] = 0.], [P[1] = 97.50000000, P[2] = 97.50000000]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 137374336), Vector[column](%id = 137697344)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 143382152), Vector[column](%id = 143472596) (17)
 

> `:=`(df, dirfield(F, 96.5, 98.5, 96.5, 98.5)); -1
 

> display(df); 1
 

Plot_2d
 

> display(df, seq(phasedia2(F, P0, 10), P0 = [`<,>`(96.5, 96.5), `<,>`(96.5, 97.5)])); 1
 

Plot_2d
 

> timeplot2(F, `<,>`(20, 30), 40); 1
 

Plot_2d
 

> timeplot2(F, `<,>`(35, 30), 40); 1
 

Plot_2d
 

>