Application Center - Maplesoft

App Preview:

Application of the Lambert W Function to the SIR Epidemic Model

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



Application of the Lambert W Function to the SIR Epidemic Model

Frank Wang, 


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, 700-721 (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, `*`(k, `*`(s(t), `*`(i(t)))).  The rate of transition from class i(t) to class r(t) is proportional to i(t), that is, `*`(l, `*`(i(t))).  With these assumptions, the differential equations describing the number of individuals in the three classes are 

diff(s(t), t) = `+`(`-`(`*`(k, `*`(s(t), `*`(i(t)))))) 

diff(i(t), t) = `+`(`*`(k, `*`(s(t), `*`(i(t)))), `-`(`*`(l, `*`(i(t))))) 

diff(r(t), t) = `*`(l, `*`(i(t))) 

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. 

> restart:

> with(plots):

> eq1 := diff(s(t),t) = -k*s(t)*i(t);

diff(s(t), t) = `+`(`-`(`*`(k, `*`(s(t), `*`(i(t)))))) (1)

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

diff(i(t), t) = `+`(`*`(k, `*`(s(t), `*`(i(t)))), `-`(`*`(l, `*`(i(t))))) (2)

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

diff(r(t), t) = `*`(l, `*`(i(t))) (3)

> i0 := 1/1000; s0 := 1-i0;


`/`(1, 1000)
`/`(999, 1000) (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);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (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 k and l 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);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (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: 

> display([ps1, ps2]);


> display([pi1, pi2]);


> display([pr1, pr2]);


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 

s[0] = `/`(`*`(l), `*`(k)) 

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, 329-359 (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. 12-22.)  

It is easy to verify that `+`(s(t), i(t), r(t)) = N, a constant.  From the first and third equations,  

s(t) = `*`(s[0], `*`(exp(`+`(`-`(`/`(`*`(k, `*`(r(t))), `*`(l))))))) 


diff(r(t), t) = `*`(l, `*`(`+`(N, `-`(r(t)), `-`(`*`(s[0], `*`(exp(`+`(`-`(`/`(`*`(k, `*`(r(t))), `*`(l))))))))))) 

Let us solve this equation 

> uf := expand(solve(1- r - s0*exp(-k*r/l) = 0, r));

`+`(`/`(`*`(l, `*`(LambertW(`+`(`-`(`/`(`*`(`/`(999, 1000), `*`(k, `*`(exp(`+`(`-`(`/`(`*`(k), `*`(l)))))))), `*`(l))))))), `*`(k)), 1) (7)

We obtain the Lambert W function. 

> uf := subs(k=xi*l, uf);

`+`(`/`(`*`(LambertW(`+`(`-`(`*`(`/`(999, 1000), `*`(xi, `*`(exp(`+`(`-`(xi)))))))))), `*`(xi)), 1) (8)

From the plot of the solution, we see a phase transition at `/`(`*`(k), `*`(l)) = 1.  (We define xi = `/`(`*`(k), `*`(l)).) 

> plot(uf, xi=0..2.5, thickness=2);


The solution is the asymptotic value of r(t).  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):

.2005682172 (9)

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

0.58982240e-2 (10)

> display([pr1, p7]);


> display([pr2, p8]);



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 s[0] = `/`(`*`(k), `*`(l)).   


In this Appendix, we examine the properties of the Lambert W function to further understand the phase-change behavior.   

> restart:

From the solution in the main text, we see that the asymptotic value of the removed class takes the form of  

`+`(`/`(`*`(LambertW(`+`(`-`(`*`(x, `*`(exp(`+`(`-`(x))))))))), `*`(x)), 1) 

Let us plot the first term: 

> f := x-> LambertW(-x*exp(-x))/x;

proc (x) options operator, arrow; `/`(`*`(LambertW(`+`(`-`(`*`(x, `*`(exp(`+`(`-`(x))))))))), `*`(x)) end proc (11)

> plot(f(x), x=0..4);


From the graph, it appears that y = -1 for `<=`(x, 1).  We plot `+`(`/`(`*`(LambertW(`+`(`-`(`*`(x, `*`(exp(`+`(`-`(x))))))))), `*`(x)), 1) below 

> plot(f(x)+1, x=0..4, thickness=2);


Below we show that `/`(`*`(LambertW(`+`(`-`(`*`(x, `*`(exp(`+`(`-`(x))))))))), `*`(x)) is indeed -1 for `<`(x, 1).  Therefore, the asymptotic value of the removed class is essentially 0 for `<`(x, 1), which explains the phase transition at x = 1.   

> simplify(f(x));

`/`(`*`(LambertW(`+`(`-`(`*`(x, `*`(exp(`+`(`-`(x))))))))), `*`(x)) (12)

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

-1 (13)

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

`/`(`*`(LambertW(`+`(`-`(`*`(x, `*`(exp(`+`(`-`(x))))))))), `*`(x)) (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 non-commercial, non-profit use only. Contact the authors for permission if you wish to use this application in for-profit activities.