Example 3: A SEIRS Spatial model

The SEIRS spatial model

Consider the population is distributed over a set of patches (regions). The individuals from one patch can move to the neighbouring ones. Consider the number of individuals in each patch is Np. The ODE represent this concern is below:

Suppose a SEIRS model in which individuals are categorised in four classes: first, all newborn individuals are assumed to be in Susceptible (S) then enter in Exposed (E) class who are infected but not yet infectious, become Infectious (I) and then change to Recovery (R ). After a period, individuals will enter to the Susceptible class again. The system of ODEs represents the SEIRS model is following:

This model can be specified through two concerns: spatial concern and SEIRS concern. It is written using the Kendrick DSL:

KModel SIR
	attribute: #(status -> S I R);
	parameters: #(beta lambda gamma mu);
	transitions: #(
		S -- lambda --> I.
		I -- gamma --> R.
		status -- mu --> Empty.
		Empty -- mu --> S.).

KModel SEIRS
	extends: 'SIR';
	parameters: #(sigma nu);
	delay: #(sigma , S -- lambda --> I  , E);
	addTransition:	#(R -- nu --> S.);
	addTransition: #(E -- mu --> Empty.).

KModel Spatial1
	attribute: #(patch -> 1  2  3  4  5);
	parameters: #(rho);
	transitions: #(
		1 -- rho --> 2.
		2 -- rho --> 3.
		3 -- rho --> 4.
		4 -- rho --> 5.
		5 -- rho --> 1.).

Composition SEIRSpatial
	model: 'Spatial1';
	model: 'SEIRS'.

Scenario SEIRSpatialPopulation
	on: 'SEIRSpatial';
	S_patch: #(900 1000 1000 1000 1000);
	I_patch: #(100 0 0 0 0).

Scenario SEIRSpatialParameters
	on: 'SEIRSpatial';
   beta_patch: #(0.75 0.5 0.5 0.5 0.5);
	N: #(patch);
	lambda: #(beta*I_patch/N);
	mu: 0.0000365;
	sigma: 0.5;
	gamma: 0.25;
	rho: 0.03;
	nu: 0.00274.
	
Simulation SEIRSSpatialRK rungeKutta
	scenarios: #(SEIRSpatialPopulation SEIRSpatialParameters);
	from: 0; 
	to: 500; 
	step: 1.
	
Visualization SEIRSSpatialViz diagram 
	for: 'SEIRSSpatialRK';
	data: #(I_patch);
	legendTitle: 'Infectious';
	xLabel: 'Time (days)';
   open.

The same model can be written also in Pharo:

spatialConcern := KEModelPart new.
spatialConcern addAttribute: #patch value: (1 to: 5) asArray.
spatialConcern addParameter: #rho.
  (1 to: 5) do: [ :i| 
        (i < 5) 
        ifTrue: [  
            spatialConcern 
                addTransitionFrom: { #patch->i } 
                to: { #patch->(i+1) } 
                probability: [ :aModel| aModel atParameter: #rho ].
        ]
        ifFalse: [  
            spatialConcern 
                addTransitionFrom: { #patch->i } 
                to: { #patch->1 } 
                probability: [ :aModel| aModel atParameter: #rho ].
        ]
    ]. 

seirsConcern := KEModelPart new.
seirsConcern attributes: {#status->#(#S #E #I #R)}.
seirsConcern addParameters: { #beta. #gamma. #mu. #sigma. #v }.
seirsConcern addParameter: #lambda value: [ :aModel | 
            (aModel atParameter: #beta) * 
            (aModel atCompartment: {#status->#I}) / (aModel atParameter: #N) ].
seirsConcern addTransitionFrom: '{#status: #S}'
        to: '{#status: #E}'
        probability: [ :aModel | 
            (aModel atParameter: #lambda) ].
seirsConcern 
        addTransitionFrom: '{#status: #E}' 
        to: '{#status: #I}' 
        probability: [ :aModel | aModel atParameter: #sigma ].
seirsConcern 
        addTransitionFrom: '{#status: #I}' 
        to: '{#status: #R}' 
        probability: [ :aModel | aModel atParameter: #gamma ].
seirsConcern 
        addTransitionFrom: '{#status: #R}' 
        to: '{#status: #S}' 
        probability: [ :aModel | aModel atParameter: #v ].
#(S E I R) do: [ :each| 
        seirsConcern 
            addTransitionFrom: {#status->each}
            to: #empty 
            probability: [ :aModel | aModel atParameter: #mu ].
        ].
seirsConcern 
        addTransitionFrom: #empty 
        to: '{#status: #S}' 
        probability: [ :aModel | aModel atParameter: #mu ].

model := spatialConcern + seirsConcern.

model atParameter: #beta assignValue: [ :aModel| |c val| 
        c := aModel currentCompartment at: #patch.
        c = 1 ifTrue: [ val := 0.75 ] ifFalse: [ val := 0.5 ].
        val 
    ].
model atParameter: #v assignValue: 0.00274.
model atParameter: #mu assignValue: 0.0000365.
model atParameter: #sigma assignValue: 0.5.
model atParameter: #gamma assignValue: 0.25.
model atParameter: #rho assignValue: 0.03. 
model atParameter: #lambda assignValue: [ :aModel| |c|
        c := aModel currentCompartment at: #patch.
        (aModel atParameter: #beta) * (aModel atCompartment: {#status->#I. #patch->c}) / (aModel atParameter: #N)
    ].

model atCompartment: { #status->#S. #patch->1 } put: 900 atOthersPut: 0.
model atCompartment: { #status->#I. #patch->1 } put: 100.
model atCompartment: { #status->#S. #patch->2 } put: 1000.
model atCompartment: { #status->#S. #patch->3 } put: 1000.
model atCompartment: { #status->#S. #patch->4 } put: 1000.
model atCompartment: { #status->#S. #patch->5 } put: 1000.
model atParameter: #N assignValue: [ :aModel| |c| 
        c := aModel currentCompartment at: #patch.
        aModel sizeOfPopulation: {c}
        ].
simulator := KESimulator new: #RungeKutta from: 0 to: 500 step: 1.
simulator executeOn: model.

chart := KEChart new.
chart
		addDataFrame: (simulator
						timeSeriesOutputsAt:
							{#status -> #I.
							#patch -> 1});
		addDataFrame: (simulator
						timeSeriesOutputsAt:
							{(#status -> #I).
							(#patch -> 2)});
		addDataFrame: (simulator
						timeSeriesOutputsAt:
							{(#status -> #I).
							(#patch -> 3)});
		addDataFrame: (simulator
						timeSeriesOutputsAt:
							{(#status -> #I).
							(#patch -> 4)});
		addDataFrame: (simulator
						timeSeriesOutputsAt:
							{(#status -> #I).
							(#patch -> 5)}).
	chart legendTitle: 'Infectious'.
	chart plot