kill(all)$

biotech : S,I incomplete model,used phase analysis,version 1.0
peter.vlasschaert@gmail.com,16/04/2020

General :compartmental modelling  (part 1)
=================================

load packages
=============

load(lrats)$

load(eigen)$

load(contrib_ode)$

Mathematical model  (S,I), incomplete
-----------------------------------------------------

symbol  : S= (s)usceptible  (= not yet infected,allowed)
                 I  = (i)nfected
S,I are : compartments (= two compartments model)

S,I  (incomplete two compartmentmodel,β,μ )
==================================
β: transmission coef.
μ: recovery  rate (= dead or be cured)

1e) compartment =S ,2e) compartment =I

info : (out) compartment  = neg sign
         (in   ) compartment =  pos sign

? :Why incomplete, because μ*I stream goes (out) system (no more part of system)

system nonlinear differential equations. (S,I)
==============================

p1a:'diff(S,t) =-β*I*S;        /*f1(x,y)=f1(S,I),ref : S compartment*/

see :simple model epidemic: dx/dt = r*x*(1-x)

p1b:'diff(I,t) = β*I*S-μ*I;   /*f2(x,y)=f2(S,I),ref; I compartemen*/

rem: S ,I = people  > 0  →   S>0,I>0
         t=time :S=S(t)    →  ic:  S(0) =S0 ,I(0)=I0
==============================================================
ic : (i)nitial (c)ondition
==============================================================
S(t) : # number susceptible   @ time t
I(t)   : # number infected         @ time t

first :why incomplete model (mathematical).
=================================

non conservative system :

p1c :p1a+p1b;

def : conservative system = rhs(p1c) = 0

Why : Incomplete model = nonconservative
=================================
→     rhs(p1c) ≠ 0

second : make system complete for S,I.
==============================

p2a:'diff(S,t) =-β*I*S+I*μ;

p2b:'diff(I,t) = β*I*S-μ*I; 

?  conservative system 
---------------------------------

p2c:p2a+p2b;

result :→  definition above, conservative 'p2c'

p2d: S+I = N;

p2d1:rhs(p2d) = lhs(p2d);

from elementary : N = constant = total population
suppose :  N =S  &  I <<<< N.

Simplification : we use nonconservative system
====================================
 → "phase plane : 2D  ,dimension  S against  I"
====================================

PHASE PLANE : ODES
===================

p3a:p1a;

p3b:p1b;

assume(S>0,I>0);

parameter : μ,β>0
------------------------

assume(β>0,μ>0);

first result :phase plane analysis
=========================

use :signum , 'command maxima'
==========================
signum(-3.0)  → -1
signum(5)       → 1

?  find sign of rhs (p3a), dS/dt
-----------------------------------------
see :assume 'above' ,why

p3d:lhs(p3a) =signum(rhs(p3a)); /*S>0 */

? find max (S,I ) plane , S>0,I>0
-------------------------------------------
everything take place : first quadrant.
euclidean space
=============
x=S
y=I

set : mathematical allowed

p3f1:p3b/p3a;

Formula : 1

p3f2:'diff(I,S) = lhs(p3f1);

p3f2a:rhs(p3f2) = lhs(p3f2);

p3f3:subst([p3f2a],p3f1),expand;

max(S(I)) ? ,  find  dI/dS = 0

p3f4:lhs(p3f3) = 0;

combination : 'p3f3,p3f4'
----------------------------------
S= smax, for above condition

p3f5:subst([p3f4,S=smax],p3f3);

peak : epidemic evolution S,I diagram

p3f6:solve(p3f5,smax);

check :max , second derivatives must be negative.
→ use p3f3

p3f7:'diff(I,S,2) = diff(rhs(p3f3),S);

→ signum ,rhs(p3f7)

p3f8:lhs(p3f7)=signum(rhs(p3f7));

final result : peak (= smax ) max in the (S,I) plane.
use : smax =threshold value
rem :not  time diagram

result : beginning 'I' very small: Infected people

p3f9a: smax= I;

Infection grows → result :can be pandemic

p3f9b:subst([p3f9a],smax > rhs(p3f6[1]));

Infection goes down →result :max limit 'every one will be recoverd'

p3f9c:subst([p3f9a],smax < rhs(p3f6[1]));

rem : I,S ,people 

what :' R0 ':
=========
use :p3f3,p3f4

def: R0 = β/μ*S

p3f10:R0 =β/μ*S;

COVID-19 : R0 >3 , normal flue:   0.5  <R0 < 2.0

What about this number:R0 (=reproduction rate)
=====================================
S : goes down by vaccination
β  : washing hands 
μ  : constant , not easy to change.(= medical care)
======================================
R0 > 1 :average grows    geometrically
R0 <  1 :average shrinks geometrically
======================================

What about this :ode , find  :S=S(I)
===========================
use:p3f3

differential equations : I(S)
------------------------------------

pg1:p3f3;

use : ode2 from 'maxima'

method 1:
=======

pg2:ode2(pg1,I,S),expand;

1e) start: value to find , %c
-------------------------------------
S(0) =S0 , I(0) = I0 ,see initial model

pg3:subst([I=I0,S=S0],pg2);

pg4:solve(pg3,%c),expand;

method 2:
========

pg5:contrib_ode(pg1,I,S)[1],expand;

check: method1,method2
====================

is(rhs(pg2) =rhs(pg5));

2e) final result  I=I(S)
----------------------------
use : pg5 & pg4

assume(S0>0,I0>0);

pg6:subst(pg4,pg5);

pg7:lhs(pg6)=logcontract(rhs(pg6));

pg8:pg7,expand;

ratprint :false;

rem : by default :ratprint = true, change option see above
reason : reduce output : decimal to fraction

intermezzo : limit computation
=======================

===================================================================

graph :log(x)
==========

wxplot2d([log(S)], [S,0,5])$

from graph : log(S) when S→0 , then log(0) = -∝  (= -infinity=minf'maxima')

find : initial : S →0 sign ? 

limit (concept,maxima)
------

pg9 :'limit(log(S),S,0)=limit(log(S),S,0);

rem :sign wrong why ? graph = -∞
                                        limit    =∞
result :sign different

S→∞

pg9a:'limit(log(S),S,∞)=limit(log(S),S,∞);

rem : wx before wxplot2d ,with this option graph inside document.

end intermezzo
========================================================================

initial value : S→0
-------------------------
see :graph

pg10:'limit('I,S,0) ='limit(rhs(pg8),S,0);

what about the sign:
use :sign, compensation,follows graph,correction=(-1)

pg10a:rhs(pg10)=(-1)*ev(rhs(pg10),nouns);

pg10b:rhs(pg10) = rhs(pg10a);

final value : S→∞
-------------------------
see :

pg11:'limit('I,S,inf) ='limit(rhs(pg8),S,inf);

pg11a:lhs(pg11)=ev(rhs(pg11),nouns);

method 1:sign:true
===========

suppose: S= 10^15  =∞
===============

pg11a1:log(10^(15)/10^2),numer;

method 2:sign
===========
use :differential equation I(S),p3f3

pg11a2:p3f3;

pg11a3:lhs(pg11a2)=limit(rhs(pg11a2),S,∞);

use ode2

pg11a4:ode2(pg11a3,I,S);

pg11a5: lhs(pg11a4) = limit(rhs(pg11a4),S,inf);

check : pg11a & pg11a5
===================

pg11a6:is(rhs(pg11a)=rhs(pg11a5));

what about this representation:
suppose :use incomplete model
================================
→ simple way , find R0 =reproduction rate
N = total population of people , can be infected (=  fixed number)
S = susceptable people
I   = infected people
=================================
see: above about conservative system : def
use :p2d,p2d1

p12a:p2d;

from elementary : N = constant = total population
suppose :  N =S  &  I <<<< N.( = true COVID19 virus)
reason : suppose I = small

1e) set : I = 0,reason N=S
===================

p12b:subst([I=0],p12a);

p12b1:rhs(p12b) =lhs(p12b);

2e) use : (S,I) incomplete model 
=========================
use : phase plane odes,'see above',for infection (=time depent),p3b

p12c:p3b;

use :p12b

p12d:subst([p12b],p12c);

3e) solve ode 12d
==============

p12e:ode2(p12d,I,t);

4e) find %c  , t=0 , I=I(t=0)=I0
======================

p12f:subst([t=0,I=I0],p12e);

p12f1:rhs(p12f)=lhs(p12f);

5e) final solution I=I(t)
================

p12g:subst([p12f1],p12e);

because that exp , may not have units
time = 1/(Nβ-μ), dimension =[s]

check solution against ode:p12d
=========================
use:p12g

p12ga:subst([p12g],p12d);

p12ga1:ev(p12ga,nouns);

p12ga2:lhs(p12ga1)=factor(rhs(p12ga1));

check: p12ga2,by hand. ( you see direct this is the same)

6e) epidemic
===========
time> 0  by exponential growth.

p12h:args(rhs(p12g));

p12h1:args(p12h[2])[2];

p2h2:subst([t = 1], p12h1)=0 ;

p2h2a:p2h2/μ,expand;

p2h3:(solve(p2h2,μ)/μ)[1];

7e) definition
==========

p2h3 :R0 = rhs(p2h3);

orginal :definition from above.
use:p12b1

p2h4:subst([p12b1],p2h3);

check : against ,def = reproduction factor= R0
use:p3f10 (=def) & p2h4

is(rhs(p3f10)=rhs(p2h4));

final result : I:f(I0,t,R0,μ)
==================

p2h5:I=I0*exp(μ*t*(R0-1));

rem :dS/dt
========

assume(R0>0);

p13:subst([p2h5],p3a);

method :1,solve ode: p13
--------------

p13a:ode2(p13,S,t);

method :2,solve ode: p13
--------------

p13b:contrib_ode(p13,S,t);

rem : p13a &p13b,same solution,see by view.

find:%c
_______
use :p13a,S(t=0)=S0

p13c:subst([t=0,S=S0],p13a);

solve :%c

p13c1:solve([p13c],%c);

final :result S
==========
use: p13b

p13c2:subst([p13c1],p13b)[1];

Model: with restriction(= S,I incomplete model, two compartment)
==================================================

with the restriction :
suppose :  N =S  &  I <<<< N.( = true COVID-19 virus)

final :combination: I & S
==================

p14da:p2h5;

p14db:p13c2;

info:S=S(t),I=I(t) 
============
S0 = S(t=0), t=time,(# susceptible & t=0)
I0   = I(t=0)               ,(# infected       &,t=0)

def :R0 =
=======
R0 : absolute values

p14dc:p3f10;

symbol  : 
=======
S  = (s)usceptible  (= not yet infected,allowed)
 I   = (i)nfected
----------------------------------------------------------------------------
S,I  (incomplete two compartmentmodel,β,μ )
==================================
β: transmission coef.
μ: recovery  rate (= dead or be cured)
-----------------------------------------------------------------------------
small : s,i
=======
s= S/N ,i=I/N , relative values  → s0= S0/N, i0=I/N
rem : N = total population (taken in account @ that time)
-----------------------------------------------------------------------------
def :r0=s*β/μ ( =relative values)
=================================================

plot : phase plane 'ODES', we need values,β,μ
=====================================
use:plotdf → construct phase plane.(p3a,p3b)
-------------------------------------------------------------------
see :plotdf , graphing for phase plane in 'Maxima'

ODES :phase plane
----------------------------

p15a:p3a;

p15b:p3b;

Where to find :β,μ (see:just above for info ,symbol)
---------------------------------------------------------------------
β,μ > 0

you tube :The SIR infectious disease model, preliminary analysis
→  https://www.youtube.com/watch?v=06wnwSEHZPY
-----------------------------------------------------------------------------------------
β  = 2
μ  =.5 ( 1/2 rocovery)
S(0) =.9 → S0
I(0)   =.1  →I0
------------------------------------------------------------------------------------------
use :p15a,p15b (relative : relation)

p15:[β=2,μ=float(1/2),S0=.9,I0=.1];

model : with values

p15a1:subst(p15,p15a);

p15b1:subst(p15,p15b);

1e) phase plane conversion model.(plotdf)
------------------------------------------------------------
use : x=I
use : y=S

p15c:[I=x,S=y];

p15c1:flatten([p15,p15c]);

x= infected,y =susceptable
------------------------------------

p15ca:subst(p15c1,p15a1);

p15cb:subst(p15c1,p15b1);

2e)  plotdf
-----------------

a) we need rhs : p15ca,15cb
------------------------------------------

p15cc:[rhs(p15ca),rhs(p15cb)];

b) plotdf
------------
smax : μ/β

p15dd:subst(p15,μ/β);

figure : global figure of the phase plane

change : xradius,yradius → smaller (see : smax)

p15cd:plotdf(p15cc,[xradius,.6],[yradius,.6])$

rem : →plateau not found : figure 3, smax
              smax = .25  , x =I = infected =.25
          →start : S (>0)→ infinity → smax → S= -infinity
=========================================
(see : above,below intermezzo )
=========================================
result (= intermezzo)
====
S →0  →  I = - infinity  (initial value)
S →∞  →  I =   infinity  (final value  )

you see change S =.25,I =.25  (max value , S,I plane)

other model : youtube
------------------------------

https://www.youtube.com/watch?v=ax0UZw8gwiE
======================================
μ =1/14 (recovery  1 to 14),β=2

time evolution 
===========
use :p14da,p14db,R0=3, t=20 days

p16:p15;

after : 4 days,time evolution

p16aa:flatten([p16,[R0=3,t=4]]);

after :20 days,time evolution

p16aa1:flatten([p16,[R0=3,t=20]]);

p16bb1:flatten([p16,[R0=3]]);

need: parameter list ,p16bb1
----------------------------------------
t = missing, draw graph :  I=I(t)

p16a:p14da;

p16aaa:subst(p16bb1,p16a);

p16aaa1:rhs(p16aaa);

wxplot2d([%], [t,0,20])$

result : 4 days (infected)

p16a1:subst(p16aa,rhs(p16a));

result : 20 days (infected)

p16a2:subst(p16aa1,rhs(p16a));

need: parameter list ,p16bb1
----------------------------------------
t=missing,draw graph : S=S(t)

p16b:p14db;

p16bbb:subst(p16bb1,p16b);

p16bbb1:rhs(p16bbb);

wxplot2d([%], [t,0,20])$

result : 4 days (susceptible)

Ħ

p16b1:subst(p16aa,rhs(p16b));

result : 20 days (susceptible)

p16b2:subst(p16aa1,rhs(p16b));

=============================================================================================