5. Several variables, SIR-model
LAST TIME

• Gave examples Fokker Planck equation of type
  
  $$\frac{\partial P(x,t)}{\partial t}=-\frac{\partial}{\partial ... ...,t)]+\frac{1}{2}\frac{\partial^2}{\partial x^2}[B(x,t) P(x,t)]$$
  
  

• A-term drift term.
• B-term represented diffusion.
• Fokker-Planck equation linear if A depends on x linearly and B is independent on x.
• Kramers equation demonstrated need to study case of more than one variable!

Today will discuss solution to linear multivariate problem.

GENERAL CASE:

$$\frac{\partial P(x_1\cdots x_r,t)}{\partial t}=-\sum_i^rc_i(t... ...m_{ij}^r B_{ij}(t)\frac{\partial^2P}{\partial x_i\partial x_j}$$

Eliminate c_i(t) term by transformation

y_i=x_i-u_i(t)

$$\Pi(y_i,t)=P(x_i-u_i(t),t)$$

Let u_i be solution to

$$\frac{du_i}{dt}=\sum_j^rA_{ij}u_j(t)+c_i(t)$$

If A_ij independent of t solution can be found by elementary methods.
In general case may need to use approximate method!
u_i(t) solution to macroscopic equation!

Find

$$\frac{\partial\Pi(y_1\cdots y_r,t)}{\partial t}= -\sum_{i,j}... ...i,j}^r B_{ij}(t)\frac{\partial^2\Pi}{\partial y_i\partial y_j}$$

Mean

$$\langle y_i\rangle =\int_{-\infty}^{\infty}y_idy_1\cdots dy_r \Pi(y_1\cdots y_r,t)$$

satisfies

$$\frac{d \langle y_i\rangle}{dt}=\sum_{i,j}^rA_{ij}\langle y_j\rangle$$

On matrix form

$$\frac{d {\bf\langle y\rangle}}{dt}={\bf A \langle y\rangle }$$

If at

$$t=0, {\bf y}={\bf y}_0$$

the solution can be expressed in terms of evolution matrix Y

$$\langle y(t)\rangle={\bf Y}(t){\bf y}_0$$

where

$$\frac{d{\bf Y}}{d t}={\bf A}(t){\bf Y}(t);\;\;{\bf Y}(0)={\bf 1}$$

Second moments

$$\frac{d\langle y_iy_j\rangle}{dt}=\sum_k^rA_{ik}\langle y_ky_j\rangle +\sum_k^rA_{jk}\langle y_iy_k\rangle+B_{ij}$$

Covariance matrix

$$C_{ij}=\langle y_iy_j\rangle-\langle y_i\rangle\langle y_j\rangle$$

satisfies same equation as second moment.
On matrix form

$$\frac{d{\bf C}}{dt}={\bf AC+CA^T+B}$$

With boundary condition

$${\bf C}(0)=0$$

solution is

$${\bf C}(t)={\bf Y}(t)\int_0^td\tau{\bf Y}(\tau){\bf B}(\tau) {\bf Y^T}(\tau)^{-1}{\bf Y^T}(t)$$

Can use this result to solve Focker-Planck equation



obtained by method of inspired guess.

Example from epidemiology
STOCHASTIC SIR MODEL

• Divide population into:
  • susceptibles S
  • infecteds I
  • immune (removed) R
• New susceptibles are introduced at the rate
  
  $$h=\Omega\gamma(1-\rho)$$
  
  

• New infecteds are introduced at a rate $\Omega\gamma\rho$. Interested in limit
  
  $$\rho<<1$$
  
  

• The susceptibles are infected and become infectious at the rate $\beta S I/\Omega$.
• The susceptibles die at the rate $\gamma S$.
• The infected individuals are removed at the rate $\lambda I$ by death or immunity.

$\Omega=$ system size parameter.

Master equation

$$\frac{\partial P(S,I,t)}{\partial t}= \sum_{S,I}[W(S,I\vert S',I')P(S',I',t)-W(S',I'\vert S,I)P(S,I,t)]$$

Where

$$W(S'I'\vert S,I)=\left\{\begin{array}{lclcl} \Omega(1-\rho)\... ... \gamma S&&S'=S-1&&I'=I\\ 0&&otherwise\\ \end{array}\right.$$

Raising and lowering operators

$$E_Sf(S,I)=f(S+1,I);\;\; E_S^{-1}f(S,I)=f(S-1,I)$$

$$E_If(S,I)=f(S,I+1);\;\;E_I^{-1}f(S,I)=f(S,I-1)$$

Express master equation in terms of raising and lowering operators

$$\frac{\partial P(S,I,t)}{\partial t}=\Omega(1-\rho)\gamma( E_... ...\+\frac{\beta}{\Omega}( E_S E_I^{-1}-1)SIP +\gamma( E_S-1)SP$$

OMEGA-EXPANSION

$$S=\Omega\phi(t)+\sqrt{\Omega}s;\;\;I=\Omega\psi(t)+ \sqrt{\Omega}i$$

Expect that $\phi$ and $\psi$ satisfy deterministic rate equations and that s and i fluctuate and satisfy a Fokker Planck equation.


P(S,I,t) probability that S susceptibles and I infecteds are present at time t.
$\Pi(s,i,t)$ probability distribution for s,i.

$$\Pi(s,i,t)=\Omega P(\Omega\phi(t)+\sqrt{\Omega} s, \Omega\psi(t)+\sqrt{\Omega}i,t)$$

Factor $\Omega$ arising from normalization condition

$$1=\sum_{S,I}P(S,I,t)=\int\int ds di \Pi(s,i,t)$$

$$=\frac{1}{\Omega} \sum_{s,i}\Pi(s,i,t)$$

We have

$$\frac{\partial P}{\partial t}=\frac{1}{\Omega}\frac{\partial\... ...{\sqrt{\Omega}}\frac{d\psi}{dt}\frac{\partial\Pi}{\partial i}$$

Next, expand

$$\frac{\partial}{\sqrt{\Omega}\partial s}+ \frac{\partial^2}{2\Omega\partial s^2}\cdots$$ E_S-1 =

$$-\frac{\partial}{\sqrt{\Omega}\partial s}+\frac{\partial^2}{2\Omega\partial s^2}\cdots$$ E_S^-1-1 =

$$\frac{\partial}{\sqrt{\Omega}\partial s}+ \frac{\partial^2}{2\Omega\partial s^2}\cdots$$ E_I-1 =

$$-\frac{\partial}{\sqrt{\Omega}\partial s}+\frac{\partial^2}{2\Omega\partial s^2}\cdots$$ E_I^-1-1 =

$$\frac{1}{\sqrt{\Omega}}\left(\frac{\partial}{\partial s}-\frac{\p... ...frac{\partial^2}{\partial s^2}-2\frac{\partial^2}{\partial s\partial i} \right)$$ E_SE_I^-1-1 =

Substitution and solving for $\frac{\partial\Pi} {\partial t}$ yields an expression of the form

$$\frac{\partial\Pi}{\partial t}=\sqrt{\Omega}\{\}+[]+O[\frac{1}{\sqrt{\Omega}}]$$

Term

$$\propto\sqrt{\Omega}$$

must vanish:

$$\{\}=\frac{\partial\Pi}{\partial s}\left(\frac{d\phi}{dt}-(1-... ...(\frac{d\psi}{dt}-\rho\gamma-\beta\psi\phi +\lambda\psi\right)$$

Get rate equations

$$\frac{d\phi}{dt}=(1-\rho)\gamma -\beta\psi\phi -\gamma\phi$$

$$\frac{d\psi}{dt}=\rho\gamma+\beta\psi\phi -\lambda\psi$$

The Fokker Planck equation is then

$$\frac{\partial\Pi}{\partial t}=+\frac{\partial}{\partial i}(\... ...ac{\partial} {\partial s}(\gamma s +\beta\psi s+\beta\phi i)\Pi$$

$$+\frac{1}{2} ([1-\rho]\gamma +\beta\psi\phi+\gamma\phi) \frac{\partial^2 \Pi}{\partial s^2}$$

$$+\frac{1}{2}(\rho\gamma+\lambda\psi+\beta\phi\psi) \frac{\pa... ...al i^2}-\beta\phi\psi\frac{\partial^2\Pi}{\partial s\partial i}$$

where solution to rate equation is substituted for $\phi$ and $\psi$.

MACROSCOPIC STEADY STATE
Put time derivatives to zero
The steady state solutions satisfy

$$0=(1-\rho)\gamma -\beta\psi\phi -\gamma\phi$$

$$0=\rho\gamma+\beta\psi\phi -\lambda\psi$$

The positive solution is stable

$$\phi_0=\frac{1}{2c}(1+c-\sqrt{(1-c)^2+4c\rho})$$

$$\psi_0=\frac{\gamma}{2c\lambda}(c-1+\sqrt{(1-c)^2+4c\rho})$$

where $c=\beta/\lambda$. If $\rho=0$ the expressions simplify to

$$\phi_0=\left\{\begin{array}{lcl} 1&for&c<1\\ \frac{1}{c}&&c>1\\ \end{array}\right.$$

$$\psi_0=\left\{\begin{array}{lcl} 0&for&c<1\\ \frac{\gamma}{\lambda}(1-\frac{1}{c})&&c>1\\ \end{array}\right.$$

We plot below $\phi$ and $\lambda\psi/\gamma$ for $\rho=.001, c=\beta/\lambda$.





When c<1 (intermittent regime) only a very small number of individuals succumb to the disease and an epidemic caused by immigration of an infectious individual will quickly die out. The population will remain disease free until the next infectious individual arrives. For c>1 the disease takes hold and there will be a finite fraction of infectious individuals in the population at all times (endemic regime).

VALIDITY OF EXPANSION

must require $\Omega\psi_0>>1,\; \Omega\phi_0>>1$.
First condition hardest to satisfy.

• If $c>1 \;\& \;\Omega>>1$ can put $\rho=0$.
• if $c<1\&\rho<<1-c$ must have
  $\gamma\rho\Omega/(\lambda(1-c))>>1$.
• if c=1 must have $\Omega\sqrt{\rho}/\lambda>>1$.

Interprete $\psi,\;\phi$ as ensemble average of many realizations. When macroscopic equations stable get Gaussian fluctuations

$$\propto\sqrt{\Omega}$$

Damped oscillation in decay of fluctuation due to drift in phase. Can show no finite-size corrections to order $\sqrt{\Omega}$.


APPROACH OF MACROSCOPIC VARIABLE TO STEADY STATE
Initially half of population infected $\rho<<1$



SUMMARY

• Can solve linear multivariate Focker-Planck equation when macroscopic equation can be solved.
• When macroscopic equations stable, fluctuations Gaussian.
• Boundedness of fluctuation provides criterion for validity of macroscopic treatemnt, in limit of large system size.
• Example SIR model from epidemiology.
• $\Omega$ exansion works when infection endemic.

Click here for Return to title page
Click here for 6. Absorbing boundaries