kill(all)$

biotech : appendix parameter μ,β (SIR model),version 1.0
peter.vlasschaert@gmail,19/04/2020

my website :biotech part
===================
see:simple model:epidemics : dx/dt=r*x*(1-x)

see: biotech : S,I incomplete model,used phase analysis,version 1.0
                        peter.vlasschaert@gmail.com,16/04/2020
====================================================
→ General :compartmental modelling  (part 1)

see biotech : S,I ,R complete model,version 1.0
                       peter.vlasschaert@gmail.com,16/04/2020
=====================================================
→ General :compartmental modelling  (part 2)

ratprint:false$

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

load(lrats)$

load(eigen)$

load(contrib_ode)$

Mathematical model  (S,I,R), complete ?
=============================================================
?  see part 1

symbol  : S= (s)usceptible  (= not yet infected,allowed)
                 I  = (i)nfected
                 R = (r)ecoverd
S,I ,R are : compartments (= three compartments model)

S,I ,R (complete three compartmentmodel,β,μ )
====================================
β: transmission coef.
μ: recovery rate (= dead or be cured)

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

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

p1a1:rhs(p1a);

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

p1b1:rhs(p1b);

p1c:'diff(R,t) = μ*I;                /*f3(x,y)=f3(S,I,R),ref; R compartement*/

p1c1:rhs(p1c);

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

info : youtube ( info :from lancet :RO)
lancet : medical yournal ( td = 6.4 days,RO = 'COVID-19' =2.68, china)
=========================================================
such time : t_halve 'radio active nuclids', radiation chemistry
                    td =doubling time ( = population doubled) =Δt =t2-t1

definition :RO (=reproduction rate = COVID-19=2.68)
=========================================
use : RO >1 ,epidemic ; RO<1,no epidemic
use :β ( = control by social distance , washing hands)

p1d:RO=β/μ*S0;

p1d1:rhs(p1d)=lhs(p1d);

find :β,μ
=======

What is doubling time (=td) : Infection
=============================

a) differential equation
-------------------------------

1e) step : dI/dt  equation
===================

p2:p1b;

2e) step: solve with 'ode2'
====================

p2a:ode2(p2,I,t);

3e) step: %c
==========
I(t=0) =I0  , find %c from initial condition

p2b:subst([t=0,I=I0,S=S0],p2a);

p2b1:rhs(p2b)=lhs(p2b);

4e) step : find general equation , I=I(t)
=============================

p2c:subst([p2b1,S=S0],p2a);

b)  How to put td in 'p2c'
----------------------------------
general info : I(t2)/I(t1) = 2 , doubling

→ It2/It1 =2, 'Maxima',formula

1e) find It1,It2
===========

p2d1:It1=subst([t=t1],rhs(p2c));

p2d2:It2=subst([t=t2],rhs(p2c));

2e) find td =t2-t1
=============
rem : log 'Maxima'  = ln 'math'

p2e: 2=rhs(p2d2)/rhs(p2d1);

p2e1:log(p2e),factor;

p2e2:subst([t2-t1=td],p2e1);

p2e3:solve(p2e2,td);

impliment :RO
---------------------
p2e3,p1d

p2e4:lratsubst([p1d1],p2e3);

c) find :μ
------------

p2e5:solve(p2e4[1],μ)[1];

use:RO=2.64,td=6.4 days, lancet, china

p2e6:[RO=2.68,td=6.4];

p2e7:subst(p2e6,p2e5);

p2e8:float(p2e7);

d) find :β
-------------
use :RO

1e) initial population :S0 
----------------------------------
N = S0+I0+R0  =approximation= S0
Belgium : S0 =population belgium =11,46 miljoen = 11460000
                  I0  = how much people infected                   (t=0)

p2f:[S0=11460000];

2e) use μ and RO
--------------------------

p2f1:p2e6;

p2f2:[p2e8];

3e) combine all information
--------------------------------------

p2f3:[p2f,p2f1,p2f2];

p2f4:flatten(p2f3);

4e) result :β
-----------------
use :definition RO,p1d

p2f5:subst([p2f4],p1d);

p2f6:solve(p2f5,β);

→ β per day

p2f7:float(p2f6[1]);

All parameters
============

p2f8:flatten([p2f4,p2f6]),float;

time evolution : differential equations for (SIR model)
============
equation : for Belgium use'p1a,p1b,p1c'

p3a:subst(p2f8,p1a);

p3a1:rhs(p3a);

p3b:subst(p2f8,p1b);

p3b1:rhs(p3b);

p3c:subst(p2f8,p1c);

p3c1:rhs(p3c);

model :rk,I0=.1,no relatives values here, 

p4aa:[p3a1,p3b1,p3c1];

inital values
=========

p4a:[S0=11460000,I0=700,R0=0];

p4a:[rhs(p4a[1]),rhs(p4a[2]),rhs(p4a[3])];

p4b:[S,I,R];

solve :wth rk ; 'Maxima'
==================
details : SIR model 'above ',before ,compartment modelling.
time :days

==============================================
syntax :
Function: rk  
    rk (ODE, var, initial, domain)  
    rk ([ODE1, …, ODEm], [v1, …, vm], [init1, …, initm], domain)

 simulation time limits
==================
[t,t_begin,t_end,day]

p4c:[t,0,200,1];

  algo
=====

sol:rk(p4aa,p4b,p4a,p4c)$

plot (SIR model)
===
time = p[1]
S,I,R = p[2],p[3],p[4]

p5a:[discrete,makelist([p[1],p[2]],p,sol)]$

p5b:[discrete,makelist([p[1],p[3]],p,sol)]$

p5c:[discrete,makelist([p[1],p[4]],p,sol)]$

wxplot2d ([p5a,p5b,p5c], [xlabel,"t=time(days)"],[ylabel,"#S,#I,#R"], [legend, "#S", "#I","#R"],grid2d)$

legend
=====

units     : μ ( 1/day),β (1/day)
#           : number of people
time      :days

simple : plot  : dI/dt ,from
===============================
→   I=I0*%e^(t*(S0*β-μ))
S0 : total population of country 
I0  : infected
see above

plot I(t),direct integration
----------------------------------

p6:I0*%e^(t*(S0*β-μ));

plot :use values from above
-------------------------------------
use:above values

p6a:flatten([p2f8,[I0=p4a[2]]]);

p6b:subst(p6a,p6);

wxplot2d(%, [t,0,200], [xlabel,"t=time(days)"] , [legend,"# infected"])$