Benchmark Model 3b

The Model M3 will be transformed to IBM model (M3b). Assuming that shedding now depends on individuals and on the time elapsed since infection, each I individual now has a transmission rate i function of the time spent in I state. The transmission rate at group scale is assumed to be the sum of individual transmission rates of group members

With Kendrick, switching from a deterministic model to a stochastic model is done just by changing the simulation entity. Kendrick supports: deterministic simulation, stochastic simulation with: Gillespie implementation, Tau-Leap Implementation (a modification of Gillespie to speed up) and IBM implementation (see Section 5, Chapter 7, page 292, Modelling infectious diseases in Humans and Animals, Keeling).

|model sirConcern ageConcern simulator db| 
model := KEModel new.

sirConcern := KEModelPart new.
sirConcern addAttribute: #status value: #(#S #I #R).
sirConcern addParameters: { #beta. #lambda. #gamma. #b. #mu }.
sirConcern 
	addTransitionFrom: '{#status: #S}' 
	to: '{#status: #I}' 
	probability: 'lambda'.
sirConcern 
	addTransitionFrom: '{#status: #I}' 
	to: '{#status: #R}' 
	probability: 'gamma'.
sirConcern
	addTransitionFrom: '{#status: #S}' 
	to: #empty 
	probability: 'mu'.
sirConcern
	addTransitionFrom: '{#status: #I}' 
	to: #empty 
	probability: 'mu'.	
sirConcern
	addTransitionFrom: '{#status: #R}' 
	to: #empty 
	probability: 'mu'.
sirConcern
	addTransitionFrom: #empty
	to: '{#status: #S}'
	probability: 'b'.

ageConcern := KEModelPart new.
ageConcern addAttribute: #age value: #(#J #A).
ageConcern addParameters: { #m }.
ageConcern 
	addTransitionFrom: '{#age: #J}' 
	to: '{#age: #A}' 
	probability: 'm'.

model integrate: sirConcern.
model integrate: ageConcern.
model population: (KEPopulation size: 100).
model addParameter: #beta_i value: [ :aModel|
	|t beta_max val|.
	t := aModel t.
	beta_max := 0.1 * (0.5 + (DhbGammaDistribution shape: 0.5 scale: 0.5) random).
	(t > 15) ifTrue: [ val := beta_max ] ifFalse: [ val := (beta_max * t / 15) ].
	val
	].
model atParameter: #beta assignValue: [ :aModel|
	|c val|
	c := aModel currentCompartment at: #age.
	c = #J ifTrue: [ val := #(1 0.01)].
	c = #A ifTrue: [ val := #(0.01 1)].
	val
	].
model atParameter: #gamma assignValue: 1/30.
model atParameter: #m assignValue: 1/50.
model atParameter: #lambda assignValue: [ :aModel|
	((aModel atParameter: #beta)*
	 (aModel atParameter: #beta_i)*
	 (aModel atCompartment: '{#status:#I}')/
	 (aModel atParameter: #N)) sum ].
model atParameter: #b assignValue: [ :aModel| 
	|c val|.
	c := aModel currentCompartment at: #age.
	c = #J ifTrue: [ val := 0.005 ].
	c = #A ifTrue: [ val := 0 ].
	val
	].
model atParameter: #mu assignValue: [ :aModel| 
	|c val|.
	c := aModel currentCompartment at: #age.
	c = #J ifTrue: [ val := 0.01 ].
	c = #A ifTrue: [ val := 0.0025 ].
	val
	].
model atCompartment: { #status->#S. #age->#A } put: 90.
model atCompartment: { #status->#I. #age->#A } put: 9.
model atCompartment: { #status->#I. #age->#J } put: 1.
simulator := KESimulator new: #IBM from: 0.0 to: 365 step: 1.
simulator executeOn: model.
db := KEDiagramBuilder new.
db data: simulator timeSeriesOutputs.
db open.