% PandoraStrain.m author: Dirk Giggenbach, see www.giggenbach.de/downloads
% start 20200402 last change 20200419
% free to copy and modify, keep the originator's name. No warranty, no official use.
% This script is only for having fun during corona-break, remaining errors are guaranteed.
% This script runs under Octave or Matlab. It simulates the pandemic Coronavirus outbrake in 2020,
% with numbers for Germany around start of April 2020. It establishes a 6-state statemachine (see "PandoraStrain.png")
% with time to state-expire and probabilities to state-change. The later will need some adjustment for "real" numbers
% since the number of "infected_not_knowing" will stay a guess until mass-tests are ongoing.
%errors: "infection_rate" and "meet_per_day" are contradictive --> change to "real_infection_rate"; see below
% One can play with different measures like testing-density or change in risk-of-contagion. Missing is seasonal effects like "lower infetion rate in summer".
% Some interesting results from this simu:
% - Total development is hardly predictable, and same start-parameters can run into favourable or really bad long-term results --> have to run this several times and average.
% This speaks for precise monitoring and modelling of the real development, to controle the "Curve" during the later development.
% - The number of "Infected_not_knowing" seems a multiple of the "Knowing_infected" in the first weeks, but these numbers get much closer after some time
% - The primary infection rate (p_H_I) was reported to be ~3 per case, but this simu indicates a number closer towards ~2.5
% - To avoid a second (even stronger) increase after the first peak (which then would overburden hospitals), secondary infection rates must be kept below 1.2 for more than a year.
% - No use saying that the picture will change completely after a vaccine is introduced.
% These results of course just apply for this simple model, which can be extended with more states or modified to better fit to reality.
%
% https://www.worldometers.info/coronavirus/country/germany/
% https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Fallzahlen.html
% https://en.wikipedia.org/wiki/2020_coronavirus_pandemic_in_Germany
%
% interesting to compare Corona numbers with last strong flu season 2017/18:
% https://edoc.rki.de/handle/176904/5739 https://www.cdc.gov/flu/about/burden/2017-2018.htm
%%% measured numbers for Germany (as of 20200406):
% tests/day: start of April = 400T/week (different sources say 350T..500T/week, seems stagnating now); of which 7% are positiv ==> means 4T positive/day
% testing: ramping up from 1 on day 1 to 31March with then 500T tests/week and 5T/d positive tests, rate staying constant at 7% positive
% (this number depends heavily on someone's decision to get tested, in reality this is not a random decision, so this might be the weakest point)
% confirmed cases (= T & M & R & D) (time): total 100132 on April 6
% confirmed in M-condition (time): - total 3936
% confirmed deaths (time): - total 1584
% confirmed recovered (time): - total 28700
% confirmend daily new cases (time): - total 28700
% 1 : healthy, might become infected
% 2 : infected and not knowing , i.e. CONTAMINOUS ; might become ill later or become tested without symptoms, or recover without symptoms after x days
% 3 : tested and infected, i.e. knowing and removed by quarantine, might become ill or dead later
% 4 : ill with symptoms i.e. detected and medical, and thus removed from beeing contaminous
% 5 : recovered
% 6 : dead
clear all; close all;
population = 8E4; % Germany-Population, in thousands
cap_lim = population/200; % limit of the medical treatment system (not discerning intensive-care from normal beds), above this number a patient dies
simu_time = 200;% 1*365; % simulating n years, assuming start of infections in Feb for Germany (first cases reported 27 Jan, but contained)
% TBD when Chinese travellers might have introduced it already earlier ...
testing_1 = 0.07/14 ; testing_2 = 1/30; % first 7% chance of positive testing per 2 weeks of conta. , then everybody is tested arbitrarily every 1 months
time_until_measures = 40; % measures at middle of March means: quarantining and more testing
time_until_relaxing = 40; % measures are relaxed after another "time_until_relaxing" causing more infections, strong testing remains
t_I_R = 14; % 2-->5 days an infected person recovers and becomes immune and non-contaminating, without treatmend and symptoms and never knowing he was infected
infection_rate = 2.0; % assumption according to authorities: each infected can infect n others in total, under normal conditions
meet_per_day_1 = 1.0; % under normal conditions, each person meets this number of others CLOSELY per day on average, and can get infected if these are from s2
meet_per_day_2 = .7; % new meets under first quarantining
meet_per_day_3 = .9; % new meets under relaxed conditions
meet_per_day = horzcat( meet_per_day_1 * ones(1, round(time_until_measures)) , ...
meet_per_day_2 * ones(1, round(time_until_relaxing )) , ...
meet_per_day_3 * ones(1, round(simu_time - (time_until_measures+time_until_relaxing) )) ) ;
p_infect = meet_per_day .* infection_rate/t_I_R ; % * s2(t)/population; % per healthy person per day , need to multiply with last factor
% error: p_infect changes here with meet_per_day*infection_rate, but one should put infection_rate=2..0,7..0,9 as stated by authorities
p_I_T = horzcat(testing_1 * ones(1, round(time_until_measures)) , testing_2 * ones(1, round(simu_time-time_until_measures)) );
% 2-->3 prob for each day that an infected gets tested ==> detected and quarantined ; then no more change in testing
t_T_R = t_I_R/2; % 3-->5 testing and quarantining happens in middle of "infected" on average
t_M_R = 21; % 4-->5 days after which an ill person recovers (unless he dies meanwhile)
p_I_M = 0.1/t_I_R; % 2-->4 likelyness of an infected to become Medical, per day
p_M_D = 0.1/t_M_R; % 4-->6 total probability that an ill person in hospital dies , per day
status = ones(population , 1); % at beginning each person is in state 1 of the six statuses
time_status = ones(population , 1); % how long has person been in current status
status(:) = 1 ; % first all persons are in status 1 = healthy on first day ...
status(1) = 2 ; time_status(1)=1; % ... except few persons are infected without knowing (only one would bear a high statistical uncertainty)
%status(2) = 2 ; time_status(2)=1;
%status(3) = 2 ; time_status(3)=1;
infected = zeros(simu_time,1);
infected(1) = 1; % each change from s1 to s2 increases the total number of infected in the population (should be equal to "population-s1(t)")
infected(2) = 3;
detected = zeros(simu_time,1); % each change s2-->s3 or s2-->s4 increases number of ever detected-infected
ever_medical = 0; % whenever a person goes to hospital
recovered = 0;
s3_sum = 0;
died = 0;
for t = 1:simu_time % counts through whole timespan given by "well_infected_p"
% disp([' Time into Simulation: ',num2str(t)]);
s1(t) = sum(status==1); % healthy
s2(t) = sum(status==2); % infected not knowing - "infected"
s3(t) = sum(status==3); % infected and tested
s4(t) = sum(status==4); % infected and ill and medical
s5(t) = sum(status==5); % recovered and immune
s6(t) = sum(status==6); % dead
allsum(t) = s1(t) + s2(t) + s3(t) + s4(t) + s5(t) + s6(t); % must always be "population"
for k = 1:population % counts through whole population
if status(k)==1 % if status is healthy (1)
if ( rand(1) < ( p_infect(t) .* (s2(t)/population) ) ) status(k)=2; infected(t)=infected(t)+1 ; time_status(k)=0; end
% likelyhood to get infected, per day, from a person in s2
elseif status(k)==2 % if status is infected & not knowing (2) these are contagious to s1
if (rand(1) < p_I_M) status(k)=4; ever_medical=ever_medical+1; detected(t)=detected(t)+1; time_status(k)=0; % Infected becomes Medical with probability per day
elseif (rand(1) < p_I_T(t)) ... % formula for tested-positive but not medical, see above
status(k)=3; s3_sum=s3_sum+1; detected(t)=detected(t)+1; time_status(k)=0; % infected become tested and thus removed from "Infected"
elseif (time_status(k) > t_I_R) status(k) = 5 ; time_status(k)=0; end
elseif status(k)==3 % if status is (3) infected and tested - they are quarantined and do not infect others anymore
if (time_status(k) > t_T_R) status(k)=5; time_status(k)=0;
end
elseif status(k)==4 % if status is (4) Medically treated (hospitalized) - they are thus quarantined and do not infect others anymore
if ( time_status(k) > t_M_R ) status(k)=5; s4(t)=s4(t)-1; recovered=recovered+1; time_status(k)=0; % Recovers after time t_M_R
elseif ( rand(1) < p_M_D ) status(k)=6; s4(t)=s4(t)-1; died=died+1; time_status(k)=0; % too long ill do die
elseif ( s4(t) > cap_lim ) status(k)=6; s4(t)=s4(t)-1; died=died+1; time_status(k)=0; % if exceeding hosüpitals' capacity, they die
end
end
time_status(k) = time_status(k) + 1; % every status is increased after one loop
end
if t/10-round(t/10)==0 disp(['on day ',num2str(t),' total ever infected (knowing and not-knowing) is ',num2str(sum(infected(:)))]); end
end
t=t;% +1; People in status_i at end of simu:
s1(t) = sum(status(:)==1);
s2(t) = sum(status(:)==2);
s3(t) = sum(status(:)==3);
s4(t) = sum(status(:)==4);
s5(t) = sum(status(:)==5);
s6(t) = sum(status(:)==6);
allsum(t) = s1(t) + s2(t) + s3(t) + s4(t) + s5(t) + s6(t); % must always be "population"
figure(100); set (gcf, "position", [gcf, gcf+200, 1000, 800]); % axis(); % numbers at the end of simulation:
ph = plot( 1 : t , s1, 'LineWidth', 2 ); hold on;
pi = plot( 1 : t , s2, 'LineWidth', 2 );
pt = plot( 1 : t , s3, 'LineWidth', 2 );
pm = plot( 1 : t , s4, 'LineWidth', 2 );
pr = plot( 1 : t , s5, 'LineWidth', 2 );
pd = plot( 1 : t , s6, 'LineWidth', 2 );
ex = plot( 1 : t , (infection_rate).^((1:t)*meet_per_day_1/t_I_R), 'LineWidth', 2 );
cap = line( [0 simu_time] , [ cap_lim cap_lim ] );
ps = plot( 1 : t , allsum, 'k-.' ); hold off; grid on;
legend( [ ph pi pt pm pr pd ex cap ] , 'Healthy' , 'Infected not knowing' , 'Tested inf. (not Med.)' , 'Medical' , 'Recovered' , 'Dead' , ...
'initial exp. growth' , ['Medical-Cap=',num2str(cap_lim)], 'Location' , 'NorthWest');
ylim([0 s5(t-1)]);
ylabel('No of persons in state 1..6'); xlabel('Time in days');
title([ 'Development of different Corona-states over time (Population=',num2str(population),')']);
tv=1:simu_time; tv_=1:simu_time-1;
real_number(tv) = infected(:);
known_number(tv) = detected(:);
real_doubling_rate(tv_) = real_number(tv_) ./ ( real_number(tv_+1) - real_number(tv_) );
known_doubling_rate(tv_) = known_number(tv_) ./ ( known_number(tv_+1) - known_number(tv_) );
figure(200); set (gcf, "position", [gcf, gcf, 800, 600]);
p_r=plot(tv_, real_number(tv_),'r') ; hold on;
p_k=plot(tv_, known_number(tv_),'b');
p_r_test=plot(tv_, population+s1(tv_),'r');
p_k_test=plot(tv_, population-s1(tv_)+s2(tv_),'b'); % see if curves match
p_d=plot(tv_, s6(tv_),'k');
hold off; grid on; title('persons infected');
legend([p_r p_k p_d], 'REAL infected', 'KNOWN infected', 'Dead', 'Location' , 'NorthWest'); axis([1 simu_time 0 max(real_number)]);
figure(300); set (gcf, "position", [gcf, gcf, 800, 600]);
p_r=plot( tv_ , real_doubling_rate(tv_) ) ; hold on; p_k=plot( tv_ , known_doubling_rate(tv_) ); grid on; hold off; ylim([0 max(real_doubling_rate)]);
legend([p_r p_k], 'REAL doubling-rate', 'KNOWN doubling-rate', 'Location' , 'NorthWest'); % axis([1 simu_time 1 max()]);
ylabel('doubling-rate in days'); title('doubling rate');
disp(' ');
disp(['Numbers on day ',num2str(simu_time),', for Population=',num2str(population),' :']);
disp(['"Healthy - never infected": ', num2str(s1(t))]);
disp(['"Infected - not knowing", on last day: ', num2str(s2(t))]);
disp(['"Tested infected" (not Med.): ', num2str(s3(t))]);
disp(['"Ever got "really infected" : ', num2str(infected(simu_time))]);
disp(['"Medically treated": ', num2str(s4(t))]);
disp(['"Medically treated", in total: ', num2str(ever_medical(t))]);
disp(['"Medically treated", Maximum: ', num2str(max(s4))]);
disp(['"Recovered" from infection: ', num2str(recovered)]);
disp(['"Officially reported Infected": ', num2str(s4(t)+s3(t)+s6(t))]);
disp(['"Died": ', num2str(s6(simu_time)),' or ',num2str(100*(died/population)),'%']);