Fit a Kendrick model to data

We give here an example on how to use an external script language like Octave (or Matlab) to fit a KENDRICK model to some data. In order to do that, we adapt the Matlab script provided by Martcheva, Maia in the book: An Introduction to Mathematical Epidemiology

Data are taken from “Influenza in a Boarding School,” British Medical Journal, 4 March 1978: https://rdrr.io/cran/outbreaks/man/influenza_england_1978_school.html we transform them in order to have 2 columns with days, number of infected

We use SIR model described page 126 of the Martcheva book:

S'(t) = - beta S(t) I(t) I'(t) = beta S(t) I(t) - alpha I(t)

Parameters of the simulation: I(3) = 25 S(3) = 763 - I(3) = 738 Population size = 763

alpha and beta parameters are pre-estimated as following: alpha = 0.3 beta = 0.0025

KENDRICK model

KendrickModel influenza_england_1978_school
	attribute: #(status -> S I R);
	parameters: #(alpha beta lambda);
	transitions: #(
		S -- lambda --> I.
		I -- alpha --> R.
	);
	lambda: #(beta*I).

Composition influenza_england_1978_school_comp
    model: 'influenza_england_1978_school'.

Scenario influenza_england_1978_school_scenario
    on: 'influenza_england_1978_school_comp';
    populationSize: 763;
    S: 738;
    I: 25;
    alpha: 0.3;
    beta: 0.0025;
    others: 0.

Simulation influenza_england_1978_school_sim rungeKutta
	scenarios: #(influenza_england_1978_school_scenario);
	from: 3.0; 
	to: 14.0; 
	step: 1;
	init;
	execute;
	exportTimeSeriesOutputsAt: '{#status:#I}'

Octave script

The script is composed of a main script and three functions. The main script run KENDRICK, fit the model to the data, print and plot the result.

% Octave script to fit a KENDRICK model to data
% Data are taken from “Influenza in a Boarding School,” British Medical Journal, 4 March 1978
% Script adapted from Marcheva book, page 126

clear all
close all

% Add Kendrick models&results in Octave path
addpath ("./Sources/Scripts/");
addpath ("./Sources/Scripts/Output/");

% loading initial data
load BSfludat.txt

% specifying higher precision 12
format long 

% define arrays with t-coordinates and y coordinates of the data
tdata = BSfludat(:,1);
qdata = BSfludat(:,2);

% Initial values of parameters alpha&beta to be fitted
alpha=0.3;
beta = 0.0025;
k = [alpha beta];

% Draw the data with initial ODE
[T Y] = ode45(@(t,y)(model_1(t,y,k)),tforward,[738.0 25.0]);
tforward = 3:0.01:14; % mesh for the solution of the differential equation
tmeasure = [1:100:1101]'; % select the points in the solution
yint = Y(tmeasure(:),2);
figure(1)
subplot(1,2,1);
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'b-');
xlabel('Time in days');
ylabel('Number of cases');
axis([3 14 0 350]);

% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
[k,fval] = fminsearch(@model_error,k);

%print final values of alpha and beta
disp(k);

%Draw the data with the final ODE
[T Y] = ode45(@(t,y)(model_1(t,y,k)),tforward,[738.0 25.0]);
yint = Y(tmeasure(:),2);
subplot(1,2,2);
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'b-');
xlabel('Time in days');
ylabel('Number of cases');
axis([3 14 0 350]);

pause(60);

% Run an iteration of the KENDRICK model with specific alpha and beta parameters
function model(v)

alpha = v(1);
beta = v(2);

%% Build KENDRICK model

model = "influenza_england_1978_school";

kendrick_model = "KendrickModel influenza_england_1978_school \
attribute: #(status -> S I R); \
parameters: #(alpha beta lambda); \
transitions: #( 		S -- lambda --> I. 		I -- alpha --> R. 	); \
lambda: #(beta*I). \
Composition influenza_england_1978_school_comp model: 'influenza_england_1978_school'.";

kendrick_model = [kendrick_model "Scenario influenza_england_1978_school_scenario on: 'influenza_england_1978_school_comp'; populationSize: 763; S: 738; I: 25; alpha:"  num2str(alpha) "; beta: " num2str(beta) "; others: 0."];

kendrick_model = [kendrick_model "Simulation influenza_england_1978_school_sim rungeKutta scenarios: #(influenza_england_1978_school_scenario); from: 3.0; to: 14.0; step: 1; init; execute; exportTimeSeriesOutputsAt: '{#status:#I}'"];

kendrick_modelname = [model ".kendrick"];

kendrick_model_filename = ["./Sources/Scripts/" kendrick_modelname];

%% Save model in file

save(kendrick_model_filename, "kendrick_model");

% Remove Octave headers for KENDRICK
system(['echo "$(tail -n +6 ' kendrick_model_filename ')" > ' kendrick_model_filename]);

%% Remove previous result of the model
system (["rm ./Sources/Scripts/Output/" model "_sim.txt"]);

%% Call KENDRICK model on the command-line
system(["./kendrick " kendrick_model_filename]);

end

% Compute the function to minimize (error between model and data)
function error = model_error(v)

load BSfludat.txt
format long
tforward = 3:0.01:14;
tmeasure = [1:100:1101]';
tdata = BSfludat(:,1);
qdata = BSfludat(:,2);

% Print values of alpha and beta parameters
disp(v);

%Run KENDRICK model
model(v);

%load file generated by KENDRICK model

load influenza_england_1978_school_sim.txt;

q = influenza_england_1978_school_sim(:,1);

% Compute the error between model and data
error = sum((q - qdata).^2);
end

% Function to solve directly the ODE in order to compare with KENDRICK result
function dy = model_1(t,y,v)

%Assign the paramets of the ODE
alpha = v(1);
beta = v(2);

dy = zeros(2,1); % assigns zeros to dy

dy(1)= -beta*y(1)*y(2); % RHS of first equation
dy(2)= beta*y(1)*y(2)-alpha*y(2); %RHS of second equation

end

In order to fit the KENDRICK model to the data, the following command should be run in the main directory of kendrick :

octave main.m

Final values of alpha and beta parameters found are : alpha = 4.650167844591154e-01 beta = 2.384159011374087e-03

The script generate the following graphs: on the left, we plot the model with the initial values of parameters and the right, you have the model plot with the final values of the parameters.