Example 2: Multi species model

A mosquito-borne disease with three host species

The SIR model with demography of a mosquito-borne disease with three host species: the primary host (mosquito) and two unspecified species, reservoir1, reservoir2. The population is partitioned using two attributes: status and species leading to 3x3 compartments. The system of ODEs representing this model is:

The multi-species concern

The transition rate matrix of this concern is zero because there is no transition from one species to another one. This concern is defined in Kendrick language as below:

multiHostConcern := KEConcern new.
multiHostConcern 
     addAttribute: #species 
     value: #(#mosquito #reservoir1 #reservoir2).

The SIR model with two species written in Kendrick language

In this script, we define directly the transitions of the model. There are 6 transitions

| model  multiHostConcern simulator db|
	multiHostConcern := KEConcern new.
	multiHostConcern 
		addAttribute: #species 
		value: #(#mosquito #reservoir1 #reservoir2).
	model := KEModel new.
	model population: (KEPopulation size: 13000).
	model attributes: {(#status -> #(#S #I #R))}.
	model addParameter: #mu value: 12.17.
	model addParameter: #gamma value: 52.
	model addParameter: #beta value: 1.
	model addParameter: #lambda value: [ :aModel| 
		(aModel atParameter: #beta) * 
		(aModel atCompartment: {#status->#I}) ].
	model
		addTransitionFrom: '{#status: #S}'
		to: '{#status: #I}'
		probability: [ :m | m atParameter: #lambda ].
	model 
		addTransitionFrom: '{#status: #I}' 
		to: '{#status: #R}' 
		probability: [ :m | m atParameter: #gamma ].
	#(#S #I #R) do: [ :each| 
		model 
			addTransitionFrom: {#status->each} 
			to: #empty 
			probability: [ :m | m atParameter: #mu ].
		].
	model 
		addTransitionFrom: #empty 
		to: '{#status: #S}' 
		probability: [ :m | m atParameter: #mu ].
	model integrate: multiHostConcern.

The multi-species concern is integrated into the SIR model. After the integration, we specify some functional transition rates of the SIR model caused by the interaction with the multi-species concern.

model 
		atParameter: #mu
		assignValue: 
		[ :aModel| |c val| 
			c := aModel currentCompartment at: #species.
			c = #mosquito ifTrue: [ val := 12.17 ].
			c = #reservoir1 ifTrue: [ val := 0.05 ].
			c = #reservoir2 ifTrue: [ val := 0.05 ].
			val
		].
	model atParameter: #N assignValue: [ :aModel| |c| 
		c := aModel currentCompartment at: #species.
		aModel sizeOfPopulation: c
		].
	model addParameter: #rho value: [ :aModel| |c val|
		c := aModel currentCompartment at: #species.
		c = #mosquito ifTrue: [ val := #(0 0.02 0.02) ].
		c = #reservoir1 ifTrue: [ val := #(0.02 0 0) ].
		c = #reservoir2 ifTrue: [ val := #(0.02 0 0) ].
		val
		].
	model atParameter: #lambda assignValue: [ :aModel|
		((aModel atParameter: #beta) * 
		(aModel atParameter: #rho) * 
		(aModel atCompartment: {#status->#I})) sum
		].

Then initialize the values for 6 compartments of the model

    model atCompartment: { #status->#S. #species->#mosquito } put: 9999.
	model atCompartment: { #status->#I. #species->#mosquito } put: 1.
	model atCompartment: { #status->#R. #species->#mosquito } put: 0.
	model atCompartment: { #status->#S. #species->#reservoir1 } put: 1000.
	model atCompartment: { #status->#I. #species->#reservoir1 } put: 0.
	model atCompartment: { #status->#R. #species->#reservoir1 } put: 0.
	model atCompartment: { #status->#S. #species->#reservoir2 } put: 2000.
	model atCompartment: { #status->#I. #species->#reservoir2 } put: 0.
	model atCompartment: { #status->#R. #species->#reservoir2 } put: 0.

The stochastic dynamics of the model using the Gillespie's direct algorithm

	simulator := KESimulator new: #Gillespie from: 0.0 to: 0.5 step: 0.0027.
	simulator executeOn: model.
	db := (KEDiagramBuilder new) data: ((simulator timeSeriesAt: '{#status: #I}') sqrt).
	db open