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 two function. The main script run KENDRICK, fit the model to the data, print and plot the result.

``octave % Octave script to fit a KENDRICK model to data % Data are provided in the Marcheva book % 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/");

load BSfludat.txt % loading initial data

format long % specifying higher precision 12

tdata = BSfludat(:,1); % define array with t-coordinates of the data

qdata = BSfludat(:,2); % define array with y-coordinates of the data

tforward = 3:0.01:14; % mesh for the solution of the differential equation

tmeasure = [1:100:1101]'; % select the points in the solution

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

[T Y] = ode45(@(t,y)(model_1(t,y,k)),tforward,[738.0 25.0]); 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);

[T Y] = ode45(@(t,y)(model_1(t,y,k)),tforward,[738.0 25.0]);

% solving the ODE with the final values of the parameters

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);

```octave
% 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

``octave % 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

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