Understand how rates are computed in Gillespie and TauLeap algorithms

[SergeStinckwich]

In the following method, the rates for the Gillespie algorithm are computed by multiplying the number of compartments by the probability. Why we are doing that?

Check if this is working for all Kendrick models (including models that are composed of 2 or 3 concerns). Add more tests/examples to check this.

doGillespieIteration: t
    | rates deltaT chosen transitions |
    rates := OrderedCollection new.
    model t: t.
    transitions := model transitions.
    transitions
        do: [ :tr | 
            | prob |
            (tr from at: #status) = #empty
                ifTrue: [ model currentCompartment: tr to ]
                ifFalse: [ model currentCompartment: tr from ].
            model nextCompartment: tr to.
            prob := (tr probability value: model) abs.
            rates add: prob * (model atCompartment: tr from) ].
    rates sum = 0
        ifTrue: [ ^ 0.0 ].
    deltaT := rand2 next ln negated / rates sum.
    chosen := self rouletteWheelSelectAmong: rates.
    (transitions at: chosen) executeOn: model times: 1.
    ^ deltaT

Same question for TauLeap algorithm:

doTauLeapIteration: t
    model t: t.
    model transitions
        do: [ :tr | 
            | rate numPoisson prob |
            "model parameters at: #scope put: tr scope.
        model parameters at: #contactingSource put: (tr from)."
            (tr from at: #status) = #empty
                ifTrue: [ model currentCompartment: tr to ]
                ifFalse: [ model currentCompartment: tr from ].
            model nextCompartment: tr to.
            prob := (tr probability value: model) abs.
            "prob isDictionary ifTrue: [ prob := prob values sum ]."
            rate := prob * (model atCompartment: tr from).
            numPoisson := (PMPoissonGenerator lambda: rate * step) next.
            tr executeOn: model times: numPoisson ]

[YvanGuifo]

The doGillespieIteration: t method returns deltaT. This method calculates the time interval for moving to a new event.

Since kendrick relies on the compartmental approach and the current version of this method has been implemented to take into account epidemiological models.

In the implementation of this method, we consider that an event can be seen as the existence of a transition. Then the rate of passage from one compartment to another is done with a probability that the event occurs.

Therefore, we consider the total rate (rates) as the product of the probability by the number of compartments

[YvanGuifo]

Can we already close this issue?