Changing initial condition for ODE

I am simulating a disease model using ODE45. What I want to achieve is that: as infection grows, I want to implement some lockdown measures that will reduce infection growth and relax the measure when infection growth as reduced to a particular value.

The system of ODE I am using is below. Which is fine.

%ODE SI model code

function [Classes] = ODE_SIR(para,ICs,maxtime)

%Run ODE using ODE45

opts = odeset('RelTol',1e-5);
[t, pop] = ode45(@diff_SIR_open_model, [1:maxtime], [ICs.S ICs.I ICs.R ICs.Sa ICs.Ia ICs.Ra], opts, para);

%Convert output to structure

Classes = struct('S',pop(:,1),'I',pop(:,2),'R',pop(:,3),'Sa',pop(:,4),'Ia',pop(:,5),'Ra',pop(:,6),'t',t);


%Diff equations

function dPop = diff_SIR_open_model(t,pop,para)

S=pop(1);
I=pop(2);
R=pop(3);
Sa=pop(4);
Ia=pop(5);
Ra=pop(6);


%The Non-adherent population
dS = - (para.beta0(1,1)*S*I + para.beta0(1,2)*S*Ia)/para.N;
dI = (para.beta0(1,1)*S*I + para.beta0(1,2)*S*Ia)/para.N -para.gamma*I ;
dR = para.gamma*I;


%The Adherent population

dSa = - (para.beta0(2,1)*Sa*I + para.beta0(2,2)*Sa*Ia)/para.Na;
dIa = (para.beta0(2,1)*Sa*I + para.beta0(2,2)*Sa*Ia)/para.Na -para.gamma*Ia;
dRa = para.gamma*Ia;


dPop = [dS; dI; dR; dSa; dIa; dRa];

end

end

However, I created a function to reduce Beta such that it impact the steady growth of the infection dynamic and reduce infection.

%% Create a function that changes rate on transmission base on NPI measure

% Input :

% Lockdown : The strenght of lockdown

% School Closure : Either school is opened or closed during the period in

% question

% Beta: The initial value of Beta when there is no control

% Output:

% Beta : The new value of Beta after control is implemented

function [new_Beta] = npi(Lockdown, School_closure,Beta)

ReductionLockdown = [1.0]; %Lockdown values

ReductionSchool = [1.0]; %School closure value

if Lockdown

ReductionLockdown = 1-Lockdown;

end

if School_closure

ReductionSchool = 1-School_closure;

end

new_Beta = Beta * (ReductionLockdown * ReductionSchool)*10;

Which is also working fine.

However when I tried to output my result, I am experiencing a sharp change in my output meanwhile I expect to see a sinuosodial kind of curve (That is the values will should not get to zero and pick back up again that is wrong). This is basically because I want to take the last value of the previous output as the next Initial condition.

How am i outputing this wrongly.

function [I, Ia,t] = npichange(para,ICs,maxtime,Beta)

%% Run model

[Classes] =ODE_SIR(para,ICs,maxtime);

for i = 1:maxtime

    if (Classes.I(end) >0) && (Classes.I(end) <=500)
        [new_Beta] = npi(0.5, 0.5,Beta);
        para.beta0 = new_Beta;
        ICs = struct('S',Classes.S(end),'I',Classes.I(end), ...
            'R',Classes.R(end),'Sa',Classes.Sa(end),'Ia',Classes.Ia(end), ...
            'Ra',Classes.Ra(end));
        [Classes] =ODE_SIR(para,ICs,maxtime);
       
    end

     if (Classes.I(end) >500) && (Classes.I(end) <=15000)
        [new_Beta] = npi(0.75, 0.75,Beta);
        para.beta0 = new_Beta;
        ICs = struct('S',Classes.S(end),'I',Classes.I(end), ...
            'R',Classes.R(end),'Sa',Classes.Sa(end),'Ia',Classes.Ia(end), ...
            'Ra',Classes.Ra(end));
        [Classes] =ODE_SIR(para,ICs,maxtime);
     end
end

t=[Classes.t; Classes.t(end)+Classes.t(2:end)]; 
        Classes.S=[Classes.S; Classes.S(2:end,:)]; Classes.Sa=[Classes.Sa; Classes.Sa(2:end,:)];
        I=[Classes.I; Classes.I(2:end,:)]; Ia=[Classes.Ia; Classes.Ia(2:end,:)];
        Classes.R=[Classes.R; Classes.R(2:end,:)]; Classes.Ra=[Classes.Ra; Classes.Ra(2:end,:)];

The values I am using to run the model

      
%Define model parameters as a structure (N.B. stick to days here for ease
%later on)
Beta = [0.5 0.1; 0.001 0.5];
para = struct('beta0',Beta,'N',0.2*62000000,'gamma',1/9, 'Na',0.8*62000000);


%Define initial conditions as a structure
ICs = struct('S',para.N-1,'I',1,'R',0,'Sa',para.Na,'Ia',0,'Ra',0);

%Define time to run model for

maxtime = 80;
[I, Ia,t] = npichange(para,ICs,maxtime,Beta);

%% Plots
figure(1)
clf
set(gca,'fontsize',16)
plot(t,I,'r',t,Ia,'k')
xlabel('Time (days)')
ylabel('Infection')
legend('Non-adherent Population', 'Adherent Population')
  

Answer 1

[SANKALP DEV]

Hi Gbeiniyi,

I understand that you are getting abrupt changes in your output.

This occurs because the values of the ‘Classes’ variable is being overwritten at each iteration of the loop. As a result, the previous values are lost, and only the values from the current iteration are retained.

To address the problem, you can store the results of each iteration separately and concatenate them at the end.

This can be done by introducing a helper function which takes the initial results and the results of each iteration, and concatenates them field by field, preserving the complete set of results.

To understand more, refer to following code:

% Helper function to concatenate the results 

function results = concatenateResults(results1,results2) 
fields = fieldnames(results1); 
for i = 1:numel(fields) 
    field = fields{i}; 
    results1.(field) = [results1.(field); results2.(field)(2:end)]; 
end 
results = results1; 
end 

%You can call the above helper function after the variable %Classes’ is being updated in the following manner: 

[Classes] = ODE_SIR(para, ICs, maxtime); 
results = concatenateResults(results, Classes);

Hope this helps you achieve the desired output!

[Gbeminiyi]

Thank you. This was very helpful