Individual-based modeling in Julia

SimJulia is a package for process-based simulations, similar to SimPy in Python, which I contributed to a little to many moons ago. Over the years, I've implemented simulation models of disease transmission to help me learn new computer languages (see my efforts in Python, C, and Eiffel on my Google Code site here). Here's my first try at a process-based version of a stochastic SIR type epidemiological model.

module SIRSim

using SimJulia
using Distributions
using DataFrames

export
    SIRPerson,
    SIRModel,
    activate,
    run,
    out

type SIRPerson
    state::Char
end

function increment(a::Array{Int64})
    push!(a,a[length(a)]+1)
end

function decrement(a::Array{Int64})
    push!(a,a[length(a)]-1)
end

function carryover(a::Array{Int64})
    push!(a,a[length(a)])
end

type SIRModel
    sim::Simulation
    parray::Array{Process}
    beta::Float64
    c::Float64
    gamma::Float64
    ta::Array{Float64}
    Sa::Array{Int64}
    Ia::Array{Int64}
    Ra::Array{Int64}
    allIndividuals::Array{SIRPerson}
end

function SIRModel(heapsize::Int64,beta::Float64,c::Float64,gamma::Float64,S::Int64,I::Int64,R::Int64)
        sim = Simulation(uint(heapsize))
        N=S+I+R
        parray = [Process(sim,string(i)) for i in 1:N]
        states = [fill!(Array(Char,S),'S'),
                  fill!(Array(Char,I),'I'),
                  fill!(Array(Char,R),'R')]
        allIndividuals=[SIRPerson(state) for state in states]
        ta=Array(Float64,0)
        push!(ta,0.0)
        Sa=Array(Int64,0)
        push!(Sa,S)
        Ia=Array(Int64,0)
        push!(Ia,I)
        Ra=Array(Int64,0)
        push!(Ra,R)
        SIRModel(sim,parray,beta,c,gamma,ta,Sa,Ia,Ra,allIndividuals)
end


function live(p::Process,individual::SIRPerson,s::SIRModel)
  while individual.state=='S'
      # Wait until next contact
      hold(p,rand(Distributions.Exponential(1/s.c)))
      # Choose random alter
      alter=individual
      while alter==individual
          N=length(s.allIndividuals)
          index=rand(Distributions.DiscreteUniform(1,N))
          alter=s.allIndividuals[index]
      end
      # If alter is infected
      if alter.state=='I'
          infect = rand(Distributions.Uniform(0,1))
          if infect < s.beta
              individual.state='I'
              push!(s.ta,now(p))
              decrement(s.Sa)
              increment(s.Ia)
              carryover(s.Ra)
          end
      end
  end
  if individual.state=='I'
      # Wait until recovery
      hold(p,rand(Distributions.Exponential(1/s.gamma)))
      individual.state='R'
      push!(s.ta,now(p))
      carryover(s.Sa)
      decrement(s.Ia)
      increment(s.Ra)
  end
end

function activate(s::SIRModel)
    [SimJulia.activate(s.parray[i],0.0,live,s.allIndividuals[i],s) for i in 1:length(s.parray)]
end

function run(s::SIRModel,tf::Float64)
    SimJulia.run(s.sim,tf)
end

function out(s::SIRModel)
    result = DataFrame()
    result["t"] = s.ta
    result["S"] = s.Sa
    result["I"] = s.Ia
    result["R"] = s.Ra
    result
end

end # module

The function, live, contains the logic behind the simulation; each susceptible individual contacts others at rate \( c \); if the contact is infected, then there is a probability \( \beta \) that the susceptible individual will become infected. After infection, individuals recover at rate \( \gamma \).

Some things to note:

• Julia doesn't permit redefining types in the main scope, so if you want to run in the REPL, types can be defined in a module, which can be reloaded.
• Julia doesn't allow subtyping of concrete types. In SimPy, one could define SIRPerson as a subclass of Process. Here, I use members of SIR model to hold instances of SimJulia.Simulation and an array of SimJulia.Process.
• I use an outer constructor method for SIRModel.
• I update S, I, and R every time step, as the output method constructs a DataFrame for the results.
• I've separated the contact rate \( c \) from the probability of infection, \( \beta \), just to clarify their different roles in transmission.
• The exclamation mark in push! denotes that the operation is in place (hence, I should really use activate! etc.).

Here is the model in action.

# Load simulation code and libraries
include("sir_sim.jl")
using SimJulia
using Gadfly

# Set parameters
# All Float64 so use decimal points
beta = 0.1
c = 1.0
gamma = 0.05

# Set initial conditions
# All Int64
S0 = 99
I0 = 1
R0 = 0
N0 = S0+I0+R0

# Initialise and run
sirmodel = SIRSim.SIRModel(N0,beta,c,gamma,S0,I0,R0);
SIRSim.activate(sirmodel);
SIRSim.run(sirmodel,1000.0);
# Collate results
result=SIRSim.out(sirmodel)

# Plot using Gadfly
p=plot(results,x="t",y="I",Geom.line)

This all seems a lot more complex than the stochastic model using Gillespie's algorithm I wrote about before. However, this approach makes it easier to build more complex models:

• The distribution of passage times (here the recovery time) can easily be changed to a distribution other than an exponential.
• Infectivity terms and contact rates can be allowed to vary in a complex way across individuals.
• Rather than assuming a well-mixed population, we can easily assume a contact network, such that individuals only contact a subset of the population.