Application of the Lambert W Function to the SIR Epidemic Model
Frank Wang, fwang@lagcc.cuny.edu
The basic framework of mathematical epidemiology is the SIR model. The letters in the acronym represent the three primary states that any member of a population can occupy with respect to a disease: susceptible, infectious, and removed. This worksheet demonstrates the use of the Lambert W function (LambertW), which was included in Maple in its early years, for the phase transition in the SIR model.
Based on the model by W. O. Kermack and A. G. McKendrick ("A contribution to the mathematical theory of epidemics" published in Proceedings of the Royal Society of London, Series A, 115, 700721 (1927)), the rate at which susceptibles become infected is proportional to the number of encounters between susceptible and infected individuals, which is proportional to the product of the two populations, . The rate of transition from class to class is proportional to , that is, . With these assumptions, the differential equations describing the number of individuals in the three classes are
It is impossible to express the solution of this system of nonlinear differential equations as explicit functions, but we can use Maple to solve them numerically.
> 
eq1 := diff(s(t),t) = k*s(t)*i(t); 

(1) 
> 
eq2 := diff(i(t),t) = k*s(t)*i(t)  l*i(t); 

(2) 
> 
eq3 := diff(r(t),t) = l*i(t); 

(3) 
> 
i0 := 1/1000; s0 := 1i0; 


(4) 
> 
soln1 := dsolve(eval({eq1, eq2, eq3, s(0)=s0, i(0)=i0, r(0)=0}, {k=1, l=0.9}), {s(t), i(t), r(t)}, numeric); 

(5) 
> 
odeplot(soln1, [t, s(t)], t=0..100); ps1 := %: 
> 
odeplot(soln1, [t, i(t)], t=0..100); pi1 := %: 
> 
odeplot(soln1, [t, r(t)], t=0..100); pr1 := %: 
We try a different set of and values:
> 
soln2 := dsolve(eval({eq1, eq2, eq3, s(0)=s0, i(0)=i0, r(0)=0}, {k=1, l=1.2}), {s(t), i(t), r(t)}, numeric); 

(6) 
> 
odeplot(soln2, [t, s(t)], t=0..100, color=green); ps2 := %: 
> 
odeplot(soln2, [t, i(t)], t=0..100, color=green); pi2 := %: 
> 
odeplot(soln2, [t, r(t)], t=0..100, color=green); pr2 := %: 
Let us plot both sets of solutions:
From the graphs, we see that in one scenario (red), there is a peak in the infected class and about 20% of the population eventually ends up in the removed class (in 1905 it meant death). In another scenario (green), both infected and removed classes are a small fraction of the total population (about 0.5% for the removed class).
Kermack and McKendrick showed that the threshold of an epidemic is at
Here we utilize Maple's LambertW function to visualize the epidemic threshold.
(For more information on the Lambert W function, see Robert M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, "On the Lambert W Function," Advances in Computational Mathematics, 5, 329359 (1996); Robert M. Corless, G. H. Gonnet, D. E. G. Hare, and D. J. Jeffrey, "Lambert's W Function in Maple," Maple Technical Newsletter 9, Spring 1993, pp. 1222.)
It is easy to verify that , a constant. From the first and third equations,
Therefore,
Let us solve this equation
> 
uf := expand(solve(1 r  s0*exp(k*r/l) = 0, r)); 

(7) 
We obtain the Lambert W function.
> 
uf := subs(k=xi*l, uf); 

(8) 
From the plot of the solution, we see a phase transition at . (We define .)
> 
plot(uf, xi=0..2.5, thickness=2); 
The solution is the asymptotic value of . Let's compare with the values from numerical solutions of differential equations.
> 
eval(uf, xi=1/.9); p7 := plot(eval(uf, xi=1/.9), t=0..100, thickness=2): 

(9) 
> 
eval(uf, xi=1/1.2); p8 := plot(eval(uf, xi=1/1.2), t=0..100, color=green, thickness=2): 

(10) 
Conclusion
In the SIR model, the asymptotic value of the removed class can be expressed as a Lambert W function. From the plot of the W function, it is clear that there is a phase transition at .
Appendix
In this Appendix, we examine the properties of the Lambert W function to further understand the phasechange behavior.
From the solution in the main text, we see that the asymptotic value of the removed class takes the form of
Let us plot the first term:
> 
f := x> LambertW(x*exp(x))/x; 

(11) 
From the graph, it appears that for . We plot below
> 
plot(f(x)+1, x=0..4, thickness=2); 
Below we show that is indeed for . Therefore, the asymptotic value of the removed class is essentially 0 for , which explains the phase transition at .

(12) 
> 
simplify(f(x)) assuming x<1; 

(13) 
> 
simplify(f(x)) assuming x>1; 

(14) 
Legal Notice: ? Maplesoft, a division of Waterloo Maple Inc. 2008. Maplesoft and Maple are trademarks of Waterloo Maple Inc. Neither Maplesoft nor the authors are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for noncommercial, nonprofit use only. Contact the authors for permission if you wish to use this application in forprofit activities.