In [35]:
using GRUtils

Viral agent simulation

Experiments with viral dynamics using an agent based simulation. Each agent is contained within an actor. There are two types of actor, a cell, which contains agents, and agents themselves. Cells are vertices in a graph and the connections, if any, between cells are the edges.

An agent may be thought of as a person or animal and a cell a physical location intimate enough that if two or more agents share a cell, they will also share the virus. Initially there are two hypotheses I would like to test with this, one is that if significant viral load can be shared between agents due to shedding, then a seeminly mild virus can suddenly wipe out a population (critical mass). Secondly that limited sampling/testing can lead to surprising jumps in reported spread. Of course this is not a remotely realistic model, it is just supposed to show some general effects observable in any complex system which resembles viral spread.

Some of the implementation is in the following external Juila script. Some structures defined there are shown below. I'm not a biologist so don't expect the names of parameters to be fully inline with the literature (putting it mildly).

If you are getting bored, just scroll down until you see graphs

In [36]:
include("virus.jl")
Out[36]:
plot_agents (generic function with 1 method)
In [37]:
@doc Agent
Out[37]:

Agent

A person or animal who may host the virus. There are many agents, each with independant parameters shown below.

  • h The current health
  • H The baselilne health
  • v The viral load
  • i The immunity to the virus
  • a Virus Antibodies
In [38]:
@doc Cell
Out[38]:

Cell

A location cosy enough for agents to share the virus

  • agents The agents in this cell
  • v The virus shedded into this cell
In [39]:
@doc World
Out[39]:

World

A container for global parameters which effect all agents and Cells.

  • R Recovery rate; how fast the agent recovers health
  • D Damage rate; effectiveness of the immune response
  • C Contagion; how fast the virus replicates within an agent regardless of health
  • CH Contagion health bias; how fast the virus replicates proportional to health
  • A Antigen; Coefficient of antiviral level
  • I Immunity; how fast immunity is developed
  • S Shedding; rate of viral shedding into the cell
In [40]:
W = World(1.025, 0.25, 1.11, 0.11, 0.9, 1.6, 0.05)
Out[40]:
World {
    R = 1.025,
    D = 0.25,
    C = 1.11,
    CH = 0.11,
    A = 0.9,
    I = 1.6,
    S = 0.05
}

Single cell sim without shedding

For each tick t of the world clock the agents are updated as follows. Note that the below function takes the current agent structure and returns a new one with updated values. This gets called by the agent actor in virus.jl. I'm sure the system of equations this represents lends itself to various mathematical analyses, but will ignore that for now and just burn some CPU cycles instead. I jiggled the parameters around until I got something that I vaguely wanted.

In [41]:
∇(w::World, a::Agent) = 
    Agent(min(w.R*a.h, a.H) - a.v - a.a, 
          a.H,
          max(a.v*(w.C + w.CH*a.h) - w.D*a.a - w.S*a.v, 0.0),
          max(a.i + w.I*a.v, 0),
          max(w.A*(a.a + a.v*a.i*a.h), 0))
Out[41]:
∇ (generic function with 1 method)

Some observations; the agent's health is damaged both by the virus and the 'antivirals'. With some viruses it is theorised that the immune response is what kills the patient in most cases, not the virus. The increase in virus can be accellerated by the agent's health, a healthy agent may have more resources for the virus or selective pressure may make the virus less aggressive when the host is already sick. To simplify things all agents develop immunity at the same rate, but healthier agent's can produce 'antivirals' quicker. Infact a healthy agent is more likely to over produce antivirals.

Frankly we could slap any number of equations together and throw some interpretation on it. It is important to understand that I'm not trying to forecast for a particular virus, instead I am just trying to observe some emergent behaviour.

In [42]:
"A message containing the starting parameters for a Single celled simulation"
struct SingleCellSim!
    W::World
    seed::Float64
end
Out[42]:
SingleCellSim!

For each tick we have some callbacks for agents and cells. The callback returns a new Agent or Cell for time t. It is called by the message handler for Tick! in virus.jl.

In [43]:
tick(::Scene, tick::Tick!{SingleCellSim!}, a::Agent) = ∇(tick.sim.W, a)
tick(::Scene, ::Tick!{SingleCellSim!}, c::Cell) = c
Out[43]:
tick (generic function with 5 methods)

Below is a message handler for an actor which creates the simulation and returns the results. In this case it just creates some agents with different baseline health values and runs a sim for 100 ticks. All the agents are put in the same cell, but they don't interact with each other because the shedding is ignored in the update rule above (although it is deducted from the viral load).

Each agent starts with a small amount of viral load and no immunity.

If you interested in how this works then invite! adds actors to the system and say/ask sends messages to actors. The actual actors are CellFrames and AgentFrames which store the Cell and Agent state at each tick.

In [44]:
function hear(s::Scene{WorldPlay}, sim::SingleCellSim!)
    wp = my(s)
    agents = [invite!(s, AgentFrames([Agent(h, h, sim.seed, 0, 0)])) for h in 0.1:0.1:1]
    wp.cells = [invite!(s, CellFrames([Cell(agents, 0)]))]
    wp.sim = sim
    wp.t = 100
    
    # Send the first tick, the Cell will tick the agents recursively and they will return tock!
    # The tocks are handled by a message handler in virus.jl which then sends out the next tick
    shout(s, wp.cells, Tick!(wp.sim, wp.t, me(s)))
end
Out[44]:
hear (generic function with 35 methods)

Start the actor system...

In [45]:
@isdefined(s) && Actors.try_say(s, stage(s), Leave!())
s, t = interact!(WorldPlay(nothing, 0, 0, [], Id(1)))
Out[45]:
(Scene{Actors.Person}, Task (runnable) @0x00007f8116c39870)

Ask the WorldPlay actor (via the Stage actor) to run the simulation by sending it a message containing the World parameters (amongst other things particular to the Actors.jl library). What we get back is an array of agent states for each tick. This is mapped into a matricies of health and virual load values and displayed. Note that the order of the columns in the matrix is semi-random because the agent data is collected asynchronously.

In [46]:
res = ask(s, stage(s), SingleCellSim!(W, 0.001), Vector);
In [47]:
res = vcat(res...) |> x -> reshape(x, 100, :)
sort!(res, dims=2; by=a -> a.H);

Plot the agent properties throughtout the simulation

In [48]:
plot_agents(res)
Out[48]:
0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Baseline health H 0 0.01 0.02 20 40 60 80 100 Time t Viral Load v 0 0.005 0.010 0.015 0.020 20 40 60 80 100 Time t Antiviral a 0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Immunity i

So the agents with a lower baseline health are killed off while the higher ones are barely effected. Agents with moderate health take much longer to recover and take longer to develop immunity. Lets try running it again, but with a much lower seed value.

In [49]:
res = vcat(ask(s, stage(s), SingleCellSim!(W, 0.0001), Vector)...) |> x -> reshape(x, 100, :)
sort!(res, dims=2; by=a -> a.H);
plot_agents(res)
Out[49]:
0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Baseline health H 0 0.005 0.010 0.015 0.020 20 40 60 80 100 Time t Viral Load v 0 0.005 0.010 0.015 0.020 20 40 60 80 100 Time t Antiviral a 0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Immunity i

This delays the onset of symptoms (visible dip in health) and reduces the amplitude of lost health. You may note that viral load goes wild for the two agents who reach a health of zero (as does immunity, but that can be ignored). I suppose this could make sense as the immune system is totally defeated. The viral load then begins to drop as the virus runs out of cells to infect and decomposition begins creating a hostile environment for the virus.

Single cell sim with shedding

Now lets try simulating spreading via shedding. I would think that in real life the amount of viral load gainded through someone elses shedded viral particulates would be too small to significantly increase your own load. Unless sharing bodily fluids in significant amounts ofcourse. However this could also model being exposed to different strains or some unknown mechanism which increases viral load or inflamation when individuals are clustered together.

Actually, as it happens, there are some papers floating about which suggest that the initial dose of virus through aerosols can have a significant effect on severity, so this is maybe not so far fetched.

In [50]:
struct SingleCellSheddingSim!
    W::World
    seed::Float64
end

struct Shed!
    v::Float64
end

Below we create some message handlers for the Shed! message above. This is a 'timed' message which is used to transfer virus from infected agents to the environment to other agents. The message is transmitted in such a way that it is guaranteed to be delivered within a given epoch (that is at time t) so that there aren't any data races.

In [51]:
function hear(s::Scene, msg::Shed!)
    frame(s).v += msg.v
end
Out[51]:
hear (generic function with 35 methods)

A new tick callback handlers are defined which sends the Shed! message with say_in which guarantees the message will be processed within a single tick. Usually messages are not guaranteed to be delivered in particular order or during a given time frame, so some extra work needs to be done to ensure that.

In [52]:
function tick(s::Scene, tick::Tick!{SingleCellSheddingSim!}, a::Agent)
    say_in(s, tick, tick.re, Shed!(a.v*tick.sim.W.S))
    
    ∇(tick.sim.W, a)
end

function tick(s::Scene, tick::Tick!, c::Cell)
    foreach(a -> say_in(s, tick, a, Shed!(c.v / 10)), c.agents)
    
    Cell(copy(c.agents), 0)
end
Out[52]:
tick (generic function with 5 methods)

For now we assume the maximum number of agents which can occupy a single cell is 10 and the proportion of transmitted virus will increase linearly as the number of agents occupying the cell increases. It is unlikely to be linear in a real 3D environment, but it may serve as a reasonable approximation.

In [53]:
function hear(s::Scene{WorldPlay}, sim::SingleCellSheddingSim!)
    wp = my(s)
    agents = [invite!(s, AgentFrames([Agent(h, h, 0, 0, 0)])) for h in 0.2:0.1:1]
    push!(agents, invite!(s, AgentFrames([Agent(0.5, 0.5, sim.seed, 0, 0)])))
    
    wp.cells = [invite!(s, CellFrames([Cell(agents, 0)]))]
    wp.sim = sim
    wp.t = 100
    
    cell = invite!(s, CellFrames([Cell(agents, 0)]))
    
    shout(s, wp.cells, Tick!(wp.sim, wp.t, me(s)))
end
Out[53]:
hear (generic function with 35 methods)
In [54]:
@isdefined(s) && Actors.try_say(s, stage(s), Leave!())
s, t = interact!(WorldPlay(nothing, 0, 0, [], Id(1)))
Out[54]:
(Scene{Actors.Person}, Task (runnable) @0x00007f80f85ff340)
In [55]:
res = vcat(ask(s, stage(s), SingleCellSheddingSim!(W, 0.001), Vector)...) |> x -> reshape(x, 100, :)
sort!(res, dims=2; by=a -> a.H);
In [56]:
plot_agents(res)
Out[56]:
0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Baseline health H 0 0.02 0.04 0.06 20 40 60 80 100 Time t Viral Load v 0 0.01 0.02 0.03 0.04 0.05 20 40 60 80 100 Time t Antiviral a 0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Immunity i

This makes things a lot worse, even the most healthy agent takes a significant dip. I suppose that is to be expected though when there are some corpses lying in the room. You may notice there is a line which begins dipping before the others and sort of waves a bit; that is the super spreader.

There is no mass extinction here, but ofcourse we could fiddle with the params until there is.

Multi-cell sim

Things can get vastly more complex moving to a mult-cell sim. We have to decide on the topology of the cells; whether they are all connected or what type of graph they form and so on. We also have to decide on agent behaviour, as they are now able to move around.

To keep things simple we will start with a fully connected graph; all cells are connected and agents can move freely between them. Imagine a spacious village with equally sized huts arranged in a circle. Each cell is a hut and the villagers teleport between them...

Also to keep things simple, agent behaviour will be random. For each tick there is some possibility that an agent will try to move to another hut picked at random, health permitting.

In [57]:
struct VillageSim!
    W::World
    "Probability any given agent will move"
    P::Float64
    "Nested array which represents cells with agents in"
    cells::Vector{Vector{Agent}}
end

struct VillageSim
    W::World
    P::Float64
    cells::Vector{Id}
end

Above you can see we now specify the cells in the simulation message as a nested vector; moving towards a declarative specification for the simulation.

In [58]:
function hear(s::Scene{WorldPlay}, simdef::VillageSim!)
    wp = my(s)
    
    
    wp.cells = [invite!(s, CellFrames([Cell([invite!(s, AgentFrames([a])) for a in agents], 0)]))
                    for agents in simdef.cells]
    wp.c = length(wp.cells)
    wp.sim = VillageSim(simdef.W, simdef.P, wp.cells)
    wp.t = 100
    
    shout(s, wp.cells, Tick!(wp.sim, wp.t, me(s)))
end
Out[58]:
hear (generic function with 35 methods)

When an agent decides to move it sends a couple of timed messages. One to the cell it wishes to move to and one to the cell it is moving from.

In [59]:
struct Ingress!
    who::Id
end

struct Egress!
    who::Id
end
In [60]:
hear(s::Scene{CellFrames}, msg::Ingress!) = push!(frame(s).agents, msg.who)
hear(s::Scene{CellFrames}, msg::Egress!) = delete!(frame(s).agents, msg.who)
Out[60]:
hear (generic function with 35 methods)

The agent decides to move in the tick callback.

In [61]:
function tick(s::Scene, tick::Tick!{VillageSim}, a::Agent)
    say_in(s, tick, tick.re, Shed!(a.v*tick.sim.W.S))
    
    if rand() > tick.sim.P && a.h > 0.3
        cell = rand(tick.sim.cells)
        
        if cell ≠ tick.re
            say_in(s, tick, cell, Ingress!(me(s)))
            say_in(s, tick, tick.re, Egress!(me(s)))
        end
    end
    
    ∇(tick.sim.W, a)
end
Out[61]:
tick (generic function with 5 methods)
In [62]:
@isdefined(s) && Actors.try_say(s, stage(s), Leave!())
s, t = interact!(WorldPlay(nothing, 0, 0, [], Id(1)))
Out[62]:
(Scene{Actors.Person}, Task (runnable) @0x00007f814db1a980)
In [63]:
sim = VillageSim!(W, 0.9, [
        [Agent(1, 1, 0.001, 0, 0), Agent(1, 1, 0, 0, 0)],
        [[Agent(h, h, 0, 0, 0), Agent(h + 0.1, h + 0.1, 0, 0, 0)] for h in 0.1:0.2:0.8]...])

res = vcat(ask(s, stage(s), sim, Vector)...) |> x -> reshape(x, 100, :)
sort!(res, dims=2; by=a -> a.H);
In [64]:
plot_agents(res)
Out[64]:
0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Baseline health H 0 0.01 0.02 0.03 0.04 20 40 60 80 100 Time t Viral Load v 0 0.01 0.02 20 40 60 80 100 Time t Antiviral a 0 0.2 0.4 0.6 0.8 1.0 20 40 60 80 100 Time t Immunity i

As you might expect this producers a more random set of curves as agents move around infecting each other.