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. In red, you have the dataset values.