SIR model with demography

\(deq\  = \ {\(Sx'\)[t]\  == \ 
          mu\ K\  - \ beta\ Sx[t]\ Ix[t]\  - \ 
            mu\ Sx[t], 
	\(Ix'\)[t]\  == \ 
          beta\ Sx[t] Ix[t]\  - \ mu\ Ix[t]\  - \ alpha\ Ix[t]\  - \ 
            gamma\ Ix[t], 
\(Rx'\)[t]\  == \ 
          gamma\ Ix[t]\  - \ mu\ Rx[t], 
Sx[0] == 
          S0, \ Ix[0] == I0, \ Rx[0] == R0}\)

Slight disease

para  =  {mu -> 0.1, beta -> 0.001,  gamma -> 0.1, 
        alpha -> 0.05,  K -> 1000}

deq /. para

deqN  =  deq /. para /. {S0 -> 1000, I0 -> 0.01,  
          R0 -> 0}

sol  =  NDSolve[deqN, {Sx[t], Ix[t], Rx[t]}, {t, 0, 100}]

Plot[Evaluate[{Sx[t], Ix[t], Rx[t]} /. sol], {t, 0, 100}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0], 
          RGBColor[0, 1, 0]}]

Plot[Evaluate[Sx[t] + Ix[t] + Rx[t] /. sol], {t, 0, 80}, 
      PlotRange -> {0, 1200}]

Mildly virulent disease

para  =  {mu -> 0.1, beta -> 0.001,  gamma -> 0.1, 
        alpha -> 0.3,  K -> 1000}

deq /. para

deqN  =  deq /. para /. {S0 -> 1000, I0 -> 0.01,  
          R0 -> 0}

sol  =  NDSolve[deqN, {Sx[t], Ix[t], Rx[t]}, {t, 0, 100}]

Plot[Evaluate[{Sx[t], Ix[t], Rx[t]} /. sol], {t, 0, 100}, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0], 
          RGBColor[0, 1, 0]}]

Plot[Evaluate[Sx[t] + Ix[t] + Rx[t] /. sol], {t, 0, 100}, 
      PlotRange -> {0, 1200}]

Severely virulent disease

para  =  {mu -> 0.1, beta -> 0.001,  gamma -> 0.1, 
        alpha -> 0.7,  K -> 1000}

deq /. para

deqN  =  deq /. para /. {S0 -> 1000, I0 -> 0.01,  
          R0 -> 0}

sol  =  NDSolve[deqN, {Sx[t], Ix[t], Rx[t]}, {t, 0, 300}]

Plot[Evaluate[{Sx[t], Ix[t], Rx[t]} /. sol], {t, 0, 300}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0], 
          RGBColor[0, 1, 0]}]

Plot[ Evaluate[Ix[t] /. sol], {t, 0, 300}, 
      PlotStyle -> RGBColor[1, 0, 0], PlotRange -> All]

Plot[Evaluate[Sx[t] + Ix[t] + Rx[t] /. sol], {t, 0, 300}, 
      PlotRange -> {0, 1200}]

eqmI  =  
      mu/beta (K - eqmS)/eqmS /. 
        eqmS -> (mu  +  alpha  +  gamma)/beta

para2  =  {mu -> 0.1, beta -> 0.001,  gamma -> 0.1,  
        K -> 1000}

eqmI /. para2

lossByDisease  =  alpha eqmI/mu

Plot[lossByDisease /. para2, {alpha, 0, 0.8}]