\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2012 (2012), No. 233, pp. 1--7.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2012 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2012/233\hfil A Lie algebra approach]
{A Lie algebra approach to susceptible-infected-susceptible epidemics}

\author[Y. Shang\hfil EJDE-2012/233\hfilneg]
{Yilun Shang}  % in alphabetical order

\address{Yilun Shang \newline
Institute for Cyber Security, University of Texas at San Antonio\\
San Antonio, Texas 78249, USA}
\email{shylmath@hotmail.com}

\thanks{Submitted August 24, 2012. Published December 21, 2012.}
\subjclass[2000]{92D30, 17B80, 60J22}
\keywords{Epidemic dynamics; Lie algebra; Riccati equation;
\hfill\break\indent susceptible-infected-susceptible}

\begin{abstract}
 The susceptible-infected-susceptible (SIS) epidemic model can be
 represented by a continuous-time Markov chain, which is governed by
 a set of deterministic differential equations (Kolmogorov forward
 equations). In this paper, a Lie algebra approach is applied to
 solve an SIS model where infection rate and recovery rate are
 time-varying. The method presented here has been used widely in
 chemical and physical sciences but not in epidemic applications due
 to insufficient symmetries.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks


\section{Introduction}

Analytical description of epidemic spreading has a long history and
can be traced back to the seminal work of Kermack and McKendrick
\cite{r1,r2}, where only three simple ordinary differential equations
are used following the mass action assumption; i.e., the rate of
increase in epidemic incidence is proportional to the product of the
number of infectious and susceptible individuals. It is also
possible to capture the propagation phenomena by mean-field theory
\cite{r3,r4,r5,r27} or generating function formalism \cite{r7,r8,r15}
especially when the host population is modeled by a network. Such
methods, however, are generally more accurate (and valid in essence)
when the population size is relatively large.

Recently, Keeling and Ross \cite{r10} proposed a time homogeneous
Markov chain model to characterize the stochastic nature of epidemic
spreading. The complete ensemble of behavior can be predicted by
$N+1$ differential equations for susceptible-infected-susceptible
(SIS) dynamics \cite{r23} by virtue of the Kolmogorov forward
equation \cite{r11}, which governs the rates of transition between
states of the disease. The solution of the system can be expressed
by the form of matrix exponentials \cite{r10,r12}. Indeed,
continuous-time Markovian models are shown to be powerful tools to
study stochastic evolutionary processes and have been widely used in
other biological and metapopulation models \cite{r14,r13,r16,r25}.

In this paper, we further investigate the SIS paradigm represented
by a time inhomogeneous Markov chain. In this model there is a fixed
population of size $N$, where $S(t)$ and $I(t)$ represent the number
of susceptibles and infectives, respectively, in the population at
time $t$, $t\ge0$, and $S(t)+I(t)=N$. No immunity is conferred upon
recovery from infection, and recovered individuals return
immediately to the susceptible state. Members of the population
transmit the infection immediately upon becoming infected. Based on
the Lie algebraic method developed in \cite{r17}, we generate a
low-dimensional Lie algebra and solve the Markovian model by
deriving matrix exponential solutions. Different from many physical
or chemical systems \cite{r19,r20}, biological or epidemic models
often lack of symmetries, which adds difficulty in finding a proper
Lie algebra. It is worth noting that Lie algebra solution of some
birth-and-death type population models is recently established by
House \cite{r18}.

The rest of the paper is organized as follows. In Section 2, we
briefly review the Lie algebraic methodology for continuous-time
Markov chains. we then apply it to an SIS epidemic model in Section
3, and conclude the paper in Section 4.


\section{Lie algebra solution of time inhomogeneous Markov chains}

In algebra theory, a Lie algebra \cite{r21} is a vector space $V$
over some field $F$ together with a bilinear map
$[\cdot,\cdot]:V\times V\rightarrow V$ called the Lie bracket, which
satisfies $[X,X]=0$ and the Jacobi identity
\begin{equation}
[X,[Y,Z]]+[Y,[Z,X]]+[Z,[X,Y]]=0,\label{e1}
\end{equation}
for all $X,Y,Z\in V$. For $X\in V$, we define an adjoint operator
$\text{ad}X$ by
\begin{equation}
(\text{ad}X)Y=[X,Y],\label{e2}
\end{equation}
for $Y\in V$. In doing so, multiple Lie brackets can be expressed in
a succinct way; e.g., $(\text{ad}X)^2Y=[X,[X,Y]]$, etc. Every
associate algebra gives rise to a Lie algebra $V$ by defining the
Lie bracket as a commutator
\begin{equation}
[X,Y]=XY-YX,\label{e3}
\end{equation}
where $X,Y\in V$. In what follows, we will focus on this Lie
product. The classical Baker-Campbell-Hausdorff formula can be
written in terms of \eqref{e2} as
\begin{equation}
e^XYe^{-X}=(e^{\text{ad}X})Y,\label{e4}
\end{equation}
where $e^X=\sum_{i=0}^{\infty}X^i/i!$.

The type of processes we consider here are continuous-time Markov
chains \cite{r10,r11}, taking values in a finite or countably infinite
state space $\mathcal{S}$. The dynamical behavior of the Markov
chain is specified by a matrix $Q(t)=(q_{ij}(t),i,j\in\mathcal{S})$,
where $q_{ij}(t)$ is the rate of transition from state $i$ to state
$j$, for $j\not=i$, and $-q_{ii}(t)=q_i(t)=\sum_{j\not=i}q_{ij}(t)$
is the total rate at which we move out of state $i$ at time $t$. In
light of the Kolmogorov forward equation (also called the ensemble
or master equation), the probability distribution of the process at
time $t$, $p(t)=(p_i(t),i\in\mathcal{S})$, is given by
\begin{equation}
\frac{{\rm d}p(t)}{{\rm d}t}=H(t)p(t),\label{e5}
\end{equation}
where $H(t)=Q(t)^T$ ($T$ means transpose), and $p(t)$ is a column
probability vector with component $p_i(t)$ representing the
probability of finding the system in state $i$ at time $t$. Making
use of the Dirac notation for vectors (kets $|\cdot\rangle$), the
probability vector can alternatively be written as
\begin{equation}
|p(t)\rangle=\sum_{i\in\mathcal{S}}P(i|t)|i\rangle,\label{e5a}
\end{equation}
where $P(i|t)$ is the probability that the Markov chain in question
taking the value of $i$ at time $t$, and $|i\rangle$ is a basis
vector, linearly independent of any other basis vector with
different value. Note that $H(t)$ in \eqref{e5} is time-dependent
implying that the process is time inhomogeneous.

The Lie algebraic method introduced in \cite{r17} requires a
decomposition of the operator $H(t)$ as
\begin{equation}
H(t)=\sum_{i=1}^ma_i(t)H_i,\label{e6}
\end{equation}
such that $a_i(t)$ are real-valued functions, and $H_i$ are linearly
independent constant operators generating a Lie algebra
$V=\text{span}\{H_1,\cdots,H_m\}$ by implementing a Lie bracket
\begin{equation}
[H_i,H_j]=H_iH_j-H_jH_i=\sum_{k=1}^m\xi_{ijk}H_k\label{e7}
\end{equation}
for some real $\xi_{ijk}$. The solution of system \eqref{e5} can be
uncoupled into a product of exponentials \cite{r17}
\begin{equation}
p(t)=e^{g_1(t)H_1}\cdots e^{g_m(t)H_m}p(0)=U(t)p(0),\label{e8}
\end{equation}
where $g_i(t)$ are real-valued functions and $g_i(0)=0$.

Substituting \eqref{e6} and \eqref{e8} into \eqref{e5}, we obtain
\begin{equation}
\begin{aligned}
\frac{{\rm d}p(t)}{{\rm d}t}
&= \sum_{i=1}^ma_i(t)H_iU(t)p(0) \\
&=\sum_{i=1}^m\dot{g}_i(t)\Big(\prod_{j=1}^{i-1}e^{g_j(t)H_j}\Big)
H_i\Big(\prod_{j=i}^{m}e^{g_j(t)H_j}\Big)p(0).
\end{aligned}\label{e9}
\end{equation}
On multiplying $U(t)^{-1}$ on both sides of \eqref{e9}, we have
\begin{equation}
\begin{aligned}
&\sum_{i=1}^ma_i(t)H_i\Big(\prod_{j=1}^me^{g_j(t)\text{ad}H_j}\Big)p(0)\\
&= \sum_{i=1}^ma_i(t)H_iU(t)p(0)U(t)^{-1} \\
&= \sum_{i=1}^m\dot{g}_i(t)\Big(\prod_{j=1}^{i-1}e^{g_j(t)H_j}\Big)
 H_i\Big(\prod_{j=i}^{m}e^{g_j(t)H_j}\Big)p(0)U(t)^{-1} \\
&= \sum_{i=1}^m\dot{g}_i(t)\Big(\prod_{j=1}^{i-1}e^{g_j(t)\text{ad}H_j}\Big)
H_i\Big(\prod_{j=1}^me^{g_j(t)\text{ad}H_j}\Big)p(0).
\end{aligned}\label{e10}
\end{equation}
Since $p(0)$ is arbitrary, we conclude that
\begin{equation}
\sum_{i=1}^ma_i(t)H_i=\sum_{i=1}^m\dot{g}_i(t)
\Big(\prod_{j=1}^{i-1}e^{g_j(t)\text{ad}H_j}\Big)H_i.\label{e12}
\end{equation}
From \eqref{e12} we derive a linear relation between $a_i(t)$ and
$\dot{g}_i(t)$ with initial values $g_i(0)=0$ (involving
$\xi_{ijk}$), as the operators $H_i$ are linearly independent.


The calculation of $p(t)$ is achieved in $O(1)$ through \eqref{e12}
rather than $O(t)$ by means of incremental direct integrations.
Therefore, the computation complexity can be dramatically reduced.
The matrix exponential form \eqref{e8} would be useful if the
derivative of the solution with respect to some model parameter is
required in subsequent calculations \cite{r17,r22}.

\section{An example: SIS epidemic spreading}

The susceptible-infected-susceptible (SIS) epidemiological model
\cite{r23} is an accurate yet simple representation of endemic
infections. It is often used as a paradigm for many sexually
transmitted infections and computer virus propagations
\cite{r2,r9,r26}. The model describes the evolution of an infection in
a fixed population, where $N$ individuals in the population are
divided into two subclasses: the susceptible pool, of size $S$, and
the infected (and infectious) class, of size $I$, with $S+I=N$.
Susceptible individuals become infected at a rate $\beta(t)$ by
contagion from infected individuals, and infected individuals, in
turn, recover (and once again become susceptible) at a rate
$\gamma(t)$.

The above description of the SIS model leads to a Markovian process
\cite{r10} whose probability vector can be written as
\begin{equation}
|p(t)\rangle=\sum_{S,I}P(S,I|t)|S,I\rangle,\label{e13}
\end{equation}
where $P(S,I|t)$ denotes the probability that there are $S$
susceptible individuals and $I$ infected ones at time $t$.
$|S,I\rangle$ is a basis vector, linearly independent of other basis
vectors with different susceptible and infected numbers. The state
space $\mathcal{S}$ consists of $N+1$ elements.

The Kolmogorov equation governing this process can be written as
\begin{equation}
\frac{{\rm d}}{{\rm d}t}|p(t)\rangle=H(t)|p(t)\rangle,\label{e15}
\end{equation}
with
\begin{equation}
H(t)=\gamma(t)(\hat{\rho}-\hat{I})+\beta(t)(\hat{\sigma}-\hat{S}),\label{e16}
\end{equation}
where
\begin{equation}
\begin{gathered}
\hat{S}|S,I\rangle=S|S,I\rangle\\
\hat{I}|S,I\rangle=I|S,I\rangle\\
\hat{\rho}|S,I\rangle=I|S+1,I-1\rangle\\
\hat{\sigma}|S,I\rangle=S|S-1,I+1\rangle
\end{gathered}\label{e17}
\end{equation}
and all these operators are linear operators (note that a similar
collection is derived for SIR model in \cite{r18}). Table 1 shows the
complete set of Lie brackets, under which the algebra
$V=\operatorname{span}\{\hat{S},\hat{I},\hat{\rho},\hat{\sigma}\}$ is
closed.

\begin{table}[ht]
\centering
\begin{tabular*}{8cm}{@{\extracolsep{\fill}}|ccccc|}
\hline
$\hat{X}$&$[\hat{X},\hat{S}]$&$[\hat{X},\hat{I}]$&$[\hat{X},\hat{\rho}]$&$[\hat{X},\hat{\sigma}]$
\\\hline
$\hat{S}$&0&0&$\hat{\rho}$&$-\hat{\sigma}$\\$\hat{I}$&0&0&$-\hat{\rho}$&$\hat{\sigma}$\\
$\hat{\rho}$&$-\hat{\rho}$&$\hat{\rho}$&0&$\hat{S}-\hat{I}$\\
$\hat{\sigma}$&$\hat{\sigma}$&$-\hat{\sigma}$&$\hat{I}$&0\\
 \hline
\end{tabular*}
\caption{Values of $[\hat{X},\hat{Y}]$ for SIS model.}
\end{table}



\begin{table}[ht]
\centering
\begin{tabular*}{10cm}{@{\extracolsep{\fill}}|ccccc|}
\hline
$\hat{X}$&$e^{g(\text{ad}\hat{X})}\hat{S}$&$e^{g(\text{ad}\hat{X})}\hat{I}$&$e^{g(\text{ad}\hat{X})}\hat{\rho}$
&$e^{g(\text{ad}\hat{X})}\hat{\sigma}$
\\\hline
$\hat{S}$&$\hat{S}$&$\hat{I}$&$e^g\hat{\rho}$&$e^{-g}\hat{\sigma}$\\
$\hat{I}$&$\hat{S}$&$\hat{I}$&$e^{-g}\hat{\rho}$&$e^{g}\hat{\sigma}$\\
$\hat{\rho}$&$\hat{S}-g\hat{\rho}$&$\hat{I}+g\hat{\rho}$&$\hat{\rho}$&$\hat{\sigma}+g\hat{S}-g\hat{I}-g^2\hat{\rho}$\\
$\hat{\sigma}$&$\hat{S}+g\hat{\sigma}$&$\hat{I}-g\hat{\sigma}$&$\hat{\rho}+g\hat{I}-\frac{g^2}{2}\hat{\sigma}$&$\hat{\sigma}$\\
 \hline
\end{tabular*}
\caption{Values of $e^{g(\text{ad}\hat{X})}\hat{Y}$ with a scalar
$g$ for SIS model.}
\end{table}


We need to look for a solution of the form
\begin{equation}
|p(t)\rangle=e^{g_1(t)\hat{S}}e^{g_2(t)\hat{I}}e^{g_3(t)
\hat{\sigma}}e^{g_4(t)\hat{\rho}}|p(0)\rangle.\label{e18}
\end{equation}
Employing \eqref{e12} and the action of exponential operators shown
in Table 2, we obtain
\begin{equation}
\begin{split}
&\gamma(t)\hat{\rho}-\gamma(t)\hat{I}+\beta(t)\hat{\sigma}-\beta(t)\hat{S}\\
&= \dot{g}_1(t)\hat{S}+\dot{g}_2(t)\hat{I}+\dot{g}_3(t)e^{g_2}e^{-g_1}
 \hat{\sigma}
+\dot{g}_4(t)\Big(e^{-g_2}e^{g_1}\hat{\rho}+g_3\hat{I}
-\frac{g_3^2}{2}e^{g_2}e^{-g_1}\hat{\sigma}\Big).
\end{split}\label{e19}
\end{equation}
Equating terms in \eqref{e19} in front of the same basis matrices
yields
\begin{equation}
\begin{gathered}
g_1(t)=-{\rm B}(t),\\
g_4(t)=\int_0^t\gamma(u)e^{g_2(u)+{\rm B}(u)}{\rm d}u,\\
\end{gathered}\label{e20}
\end{equation}
where ${\rm B}(t)=\int_0^t\beta(u){\rm d}u$; $g_2(t)$ and $g_3(t)$
are determined by the initial value problem
\begin{equation}
\begin{gathered}
\dot{g}_2(t)=-\gamma(t)-g_3(t)\gamma(t)e^{g_2(t)+{\rm B}(t)},\\
\dot{g}_3(t)=\beta(t)e^{-{\rm B}(t)-g_2(t)}
+\frac{\gamma(t)}{2}e^{g_2(t)+{\rm B}(t)}g_3(t)^2,\\
g_2(0)=g_3(0)=0.
\end{gathered} \label{e21}
\end{equation}
The function $g_3$ satisfies a Riccati equation, which may be solved
by standard reduction techniques or numerical integration; see e.g.
\cite{r24}. Let $|\mathcal{I}(t)\rangle=\sum_{S,I}I|S,I\rangle$, and
then $I(t)=\langle\mathcal{I}(t)|p(t)\rangle$ is the number of
infected individuals in the population at time $t$. In Fig. 1 we
illustrate $I(t)$ with respect to different choices of $\beta(t)$
and $\gamma(t)$ in a population of size $N=100$.


\begin{figure}[htb]
\includegraphics[width=0.7\textwidth]{fig1}
\caption{Dynamics of the SIS model with $N=100$ and
$|p(0)\rangle=|0.9\cdot N,0.1\cdot N\rangle$.
$I(t)=\langle\mathcal{I}(t)|p(t)\rangle$ are plotted with respect to
$\beta(t)=1,\gamma(t)=0.5$ (circles), $\beta(t)=0.5,\gamma(t)=1$
(squares), and $\beta(t)=e^t,\gamma(t)=1$ (triangles), where
$|p(t)\rangle$ are obtained from \eqref{e18}.}
\end{figure}

When $\gamma=0$, the model reduces to a simple susceptible-infected
(SI) epidemics, where individuals, once infected, are infected (and
infectious) forever. In this case, the solution of \eqref{e19} can be
obtained as
\begin{equation}
\begin{gathered}
g_1(t)=-{\rm B}(t),\\
g_3(t)=\int_0^te^{-{\rm B}(u)}\beta(u){\rm d}u,\\
g_2(t)=g_4(t)=0,
\end{gathered}\label{e23}
\end{equation}
where ${\rm B}(t)=\int_0^t\beta(u){\rm d}u$. Note that this can be
derived similarly from the SIR model addressed in \cite{r18}. The
consistency confirms that our result is valid.

\subsection*{Conclusion}

In this paper, we showed that it is possible to solve
susceptible-infected-susceptible (SIS) model via a Lie algebra
methodology. Lie algebra solution of differential equations has
found host of useful applications in physical systems, where wealthy
symmetries exist. Due to insufficient symmetry, this method is not
widely used in biological or social systems. For future work, more
complex and realistic epidemiological mechanisms, such as
susceptible-exposed-infectious-recovered (SEIR) model, are worthy of
further research.




\begin{thebibliography}{00}
\bibitem{r14} D. Alonso, A. McKane; 
\emph{Extinction dynamics in
mainland-island metapopulations: an N-patch stochastic model}, Bull.
Math. Biol. 64(2002) 913--958.

\bibitem{r19} A. Alvermann, H. Fehske; 
\emph{High-order commutator-free exponential time-propagation of driven
quantum systems}, J. Comput. Phys. 230(2011) 5930--5956.

\bibitem{r2} R. M. Anderson, R. M. May; 
\textit{Infectious Diseases of
Humans}, Oxford University Press, Oxford, 1991.

\bibitem{r13}A. D. Barbour, M. J. Luczak; 
\emph{A law of large numbers approximation
for Markov population processes with countably many types}, Probab.
Theory Relat. Fields 153(2012) 727--757.

\bibitem{r4}M. Barth\'el\'emy, A. Barrat, R. Pastor-Satorras, A. Vespignani; 
\emph{Dynamical patterns of epidemic outbreaks in complex
heterogeneous networks}, J. Theor. Biol. 235(2005) 275--288.

\bibitem{r20}S. Blanes, F. Casas, J. A. Oteo, J. Ros; 
\emph{The Magnus expansion and some of its
applications}, Phys. Rep. 470(2009) 151--238.

\bibitem{r5}S. G\'omez, J. G\'omez-Garde{\~n}es, Y. Moreno, A. Arenas; 
\emph{Nonperturbative heterogeneous mean-field approach to
epidemic spreading in complex networks}, Phys. Rev. E 84(2011)
036105.

\bibitem{r18}T. House; 
\emph{Lie algebra solution of population models based
on time-inhomogeneous Markov chains}, J. Appl. Probab. 49(2012)
472--481.

\bibitem{r21} J. E. Humphreys; 
\emph{Introduction to Lie Algebras and Representation
Theory}, Springer, New York, 1972.

\bibitem{r10} M. J. Keeling, J. V. Ross; 
\emph{On methods for studing stochastic disease dynamics}, 
J. R. Soc. Interface 5(2008) 171--181.

\bibitem{r1} W. O. Kermack, A. G. McKendrick; 
\emph{A contribution to the
mathematical theory of epidemics}, Proc. R. Soc. Lond. A 115(1927)
700--721.

\bibitem{r7} M. E. J. Newman; 
\emph{Spread of epidemic disease on
networks}, Phys. Rev. E 66(2002) 016128.
\bibitem{r11}J. R. Norris; \textit{Markov Chains}, Cambridge University
Press, New York, 1997.

\bibitem{r8} M. Marder; 
\emph{Dynamics of epidemics on random networks}, Phys. Rev. E 75(2007) 066103.

\bibitem{r16} O. Ovaskainen; 
\emph{The quasistationary distribution of the stochastic logistic
model}, J. Appl. Probab. 38(2001) 898--907.

\bibitem{r3} R. Pastor-Satorras, A. Vespignani; 
\emph{Epidemic spreading in scale-free networks}, Phys. Rev. Lett. 86(2001)
3200--3203.

\bibitem{r24} A. D. Polyanin, V. F. Zaitsev; 
\emph{Handbook of Exact Solutions for
Ordinary Differential Equations}, 2nd Edition, CRC Press, Boca
Raton, 2003.

\bibitem{r12}J. V. Ross; 
\emph{Computationally exact methods for stochastic
periodic dynamics spatiotemporal dispersal and temporally forced
transmission}, J. Theor. Biol. 262(2010) 14--22.

\bibitem{r25}Y. Shang; 
\emph{Likelihood estimation for stochastic epidemics
with heterogeneous mixing populations}, Int. J. Comput. Math. Sci.
6(2012) 34--38.

\bibitem{r15}Y. Shang; 
\emph{Distribution dynamics for SIS model on random
networks}, J. Biol. Syst. 20(2012) 213--220.

\bibitem{r9}Y. Shang; 
\emph{Optimal attack strategies in a dynamic botnet defense
model}, Appl. Math. Inf. Sci. 6(2012) 29--33.

\bibitem{r26}Y. Shang; 
\emph{Optimal control strategies for virus spreading in
inhomogeneous epidemic dynamics}, Canad. Math. Bull.
doi:10.4153/CMB-2012-007-2.

\bibitem{r27}Y. Shang; 
\emph{Mixed SI(R) epidemic dynamics in random graphs with
general degree distributions}, Appl. Math. Comput.
doi:10.1016/j.amc.2012.11.026.

\bibitem{r17} J. Wei, E. Norman; 
\emph{Lie algebra solution of linear
differential equations}, J. Math. Phys. 4(1963) 575--581.

\bibitem{r23} G. H. Weiss, M. Dishon; 
\emph{On the asymptotic behavior of the
stochastic and deterministic models of an epidemic}, Math. Biosci.
11(1971) 261--265.

\bibitem{r22} R. M. Wilcox; \emph{Exponential operators and parameter
differentiation in quantum physics}, J. Math. Phys. 8(1967)
962--982.

\end{thebibliography}
\end{document}