[In this module, students will become familiar with Markov Chain Monte Carlo compartmental modelling by sampling from the Exponential distribution to simulate the time between state transitions. Several examples will be given, including modelling of Aedes aegypti populations, the vector for Zika virus.]
In this past module, we discussed the Exponential and Poisson distributions. In particular we noted the duality of these two distributions; the Exponential distribution describes the time between events in a Poisson process.
In this past module, we discussed how to perform simple Markov Chain Monte Carlo stochastic modelling of a compartmental model, by sampling at each small time step the change in compartments from the Poisson distribution.
Because of the duality of the Exponential and Poisson distributions we can just as easily perform MCMC stochastic compartmental modelling randomly Exponentially sampling the time between state transitions.
In order to do this, we make use of this special property of the Exponential distribution: If X~Exp(gamma) and Y~Exp(rho), then min(X,Y) is distributed as Exp(lambda) with lambda=gamma+rho
The proof of this property can be found here.
Thus, for example, given a model with two possible state transitions, with expected rate gamma and rho, respectively, we expect the rate for at least one of the transitions to take place to be gamma+rho. Thus, we would randomly sample X~Exp(gamma) and Y~Exp(rho), and whichever is smallest is the state transition we allow to take place, and we increment the time by delta_t = min(X,Y).
Take, for example an SIR model: 
The list of state transitions is:
Thus, infection is expected to occur at rate beta*S*I/N, and recovery is expected to occur at rate gamma*I. Thus, we would randomly sample X~Exp(beta*S*I/N) and Y~Exp(gamma*I); if X is less than Y then an infection takes place, and we increment time by delta_t=X. If Y is less than X, then a recovery takes place, and delta_t=Y.
The R script sir_markov_chain_iterated_exponential.R performs an MCMC stochastic simulation of this system, with options to sample from the Exponential distribution (lexponential=1 in the script), or the Poisson distribution (lexponential=0). The script produces the plot:
Try running the script again with lexponential=0 to see how the plot changes.
The Susceptible, Exposed, Infected, Recovered model is described by the set of differential equations:
The table of possible transitions is
In this case, at each step of the simulation we randomly sample X~Exp(beta*S*I/N), Y~Exp(kappa*E), and W~Exp(gamma*I), and find the minimum, and select that as the transition that occurred.
The script seir_markov_chain_iterated_exponential.R performs 50 realisations of a MCMC stochastic SEIR model, sampling from the Exponential probability distribution to simulate the time between state transitions. The script produces the following plot: 
Dengue virus is a mosquito borne arbovirus, primarily spread between humans by Aedes aegypti, a tropical/semi-tropical mosquito. Because of the wide range of the mosquito across the world, approximately one third to half of the world’s population is at risk for DENV.
DENV typically causes a flu-like illness, but in some cases can cause haemorrhagic fever and death. There are four serotypes of the virus, and infection with one will confer immunity to it, but not the other three. Thus it is possible to get DENV several times. There are currently no vaccines or specific treatment for DENV.
DENV does not spread directly from human-to-human. Instead, a mosquito bites an infected person, and after a period of time (known as the extrinsic incubation period) the virus replicates in the mosquito’s gut, such that when it takes another blood meal it infects another person.
Chikungunya virus, CHIKV, is also spread Aedes aegypti mosquitoes. Unlike DENV it only has one known serotype, so recovery from CHIKV confers lasting immunity. However, new susceptible babies are always being born into the population, thus periodic outbreaks can occur when enough susceptibles build up, and an infected person is introduced to the population. CHIKV rarely kills, but can cause crippling sequelae, including severe arthritis. There is currently no vaccine or specific treatment for CHIKV. In 2014/15 there was an explosive outbreak of CHIKV in the Americas; the entire population was susceptible at the beginning of the outbreak.
Zika virus is also spread by Aedes. Like CHIKV, there is only one serotype, so recovery results in lasting immunity. There is currently no vaccine or specific treatment for Zika. Compared to DENV and CHIKV, Zika results in the most mild infection, and almost never kills. In fact, it has been estimated that 80% of cases are asymptomatic. However, there is recent evidence from Brazil and other countries that strongly suggests that infection with Zika virus early in pregnancy can cause significant birth defects in the fetus. Thus, the most telling evidence of a large Zika outbreak is not the burden on the healthcare system of those who are actively infected, but the babies born with defects 6 to 9 months later.
It is for this reason that the WHO has declared a global health emergency due to Zika at the beginning of February, 2016. Particularly worrisome is the fact that the 2016 summer Olympics are to be held in Brazil, meaning that there will likely be a huge influx of largely susceptible people to the region, and highly elevated probabilities of worldwide dispersal of the virus. The only bright side is that the peak season for Aedes aegypti in Brazil is January to May...
It is worth noting here that Zika comes under the category of “neglected tropical diseases”, which are diseases primarily affecting low-income countries, and get very little research funding. Zika has been so neglected, that very little study has been done of its epidemiology, including its incubation and infectious period, and the extrinsic incubation period in mosquitoes.
Currently, the only known means to control the spread of these viruses is to either prevent the mosquitoes from biting by wearing long sleeves or repellents (bed nets are not much use because Aedes primarily bite during the day), or to control the mosquito populations through the use of larvacides, insecticides, or removing breeding containers. Almost complete eradication of Aedes aegypti in the Americas was achieved in the 1960’s with the latter methods by a man named Fred Soper. Between the 1930’s to the 1960’s using DDT, and exhaustive house breeding-site searches, he managed to virtually eliminate Aedes aegypti throughout the Americas, and with it Dengue and Yellow Fever.
Aedes aegypti are “domestic” mosquitoes. They had evolved along with humans, and prefer to breed in the nice clean water we store in tanks in our houses rather than in water in rivers or ponds. We are also their favourite snack. A tank of water can only hold so many larvae, however, because they are cannabilistic. Perhaps largely due to this, the number of Aedes aegypti per household is usually quite modest… usually just a few per person. Note that only female mosquitoes bite… they need a blood meal in order to lay eggs.
Mosquitoes tend not to travel very far during their lifetimes. For instance, an analysis of female mosquitoes who took blood meals at animals at a zoo found that they travelled on average about 100 metres.
This study examined the dispersal of Aedes aegypti males, and also found that their average dispersal was around 100 m.
It is believed that the spread of diseases like DENV, CHIKV and Zika is less due to dispersion of mosquitoes, and more due to the movement of humans between patches of mosquitoes.
As an aside, and related to dispersal distances; why do you think Aedes are more prevalent in cities?
A novel means of Aedes control involves Wolbachia bacteria. Wolbachia are bacteria that live within insect cells and are passed from one generation to the next through the insect’s eggs. Wolbachia is present in up to 60% of all the different species of insects around us including some mosquitoes that bite people, but not Aedes aegypti – the primary mosquito species involved in the transmission of DENV, CHIKV, and Zika.
We can produce Wolbachia-infected Aedes aegypti mosquitoes artificially in the laboratory (let’s refer to them as Wolbachia mosquitoes). It turns out that Wolbachia mosquitoes have some interesting properties…. infection with Wolbachia impedes the the ability of DENV (and presumably CHIKV, and Zika) to replicate in the tissues of the mosquitoes.
Also, when a male mosquito infected with Wolbachia mates with a wild female mosquito, the wild female will not produce viable offspring. When Wolbachia female mosquitoes mate with either Wolbachia males or wild males, they produce Wolbachia mosquitoes.
There is a hitch though; infection with Wolbachia appears to make Aedes mosquitoes live shorter lives. Also, it appears that their larvae aren’t quite as cannabilistic as wild larvae. Also, their eggs are more likely to die if dessicated. This reduced fitness is almost certainly why we don’t see Wolbachia mosquitoes naturally in the wild.
However, the two very nice properties of Wolbachia mosquitoes (suppression of DENV and reduction in viable progeny of wild females) are so hopeful that the Bill and Melinda Gates Foundation has invested large amounts of money to examine the potential further. The effect on wild populations of releases of all male Wolbachia mosquitoes have been studied, as well as the effect of releases of both male and female Wolbachia mosquitoes. Note that the problem with releasing females is that even though they aren’t good transmitters of DENV, they still bite, and that annoys the populace.
Researchers in Australia did a test project where they did 10 releases, over 10 weeks, of male and female Wolbachia mosquitoes in the areas of about 1200 houses in Cairns, Australia. Each week, approximately 20 to 25 mosquitoes per house were released. Before the releases began, the researchers went door to door and got people to empty all sources of water in their yards.
As an aside here, Wolbachia isn’t always a good thing when it comes to human diseases; the disease filariasis is spread by a parastic worm that depend on Wolbachia to grow. Controlling that disease requires reducing Wolbachia in the worms.
If we wish to use mathematical models to predict what will happen when Wolbachia mosquitoes are released, we need to take into account the inherently patchy nature of mosquito populations. There may be, for instance, 50,000 houses in a city, with on average 10 wild Aedes per house. But that doesn’t mean that every house has Aedes, and it certainly doesn’t mean that all 500,000 mosquitoes mix together and mate. They are much more likely to mate with mosquitoes in near geographic proximity. And things like roads, walls around houses, or industrial areas, can form boundaries of their patch that are not impossible to cross, but definitely represent barriers that impede their dispersal.
Because there are so few mosquitoes per house, and geographic patches like a city block or clusters of houses aren’t that large, stochastic models are appropriate. In the R file wolbach.R, I provide an example of a stochastic MCMC model for populations of mixed wild and Wolbachia mosquitoes. The model I use has to take into account that there may be unequal numbers of males and females, particularly in the Wolbachia compartment, and that Wolbachia males and wild females cannot produce viable offspring. The model I use in the script is (note the X’s are the wild mosquitoes, and the W’s are the Wolbachia mosquitoes):
Where A_Xf is the number of babies in the X compartment due to pairings of Xf, A_Xm is the number of babies in the X compartment due to pairings of Xm, and so on….
And where B_Xf is the number of babies in the W compartment due to pairings of Xf, B_Xm is the number of babies in the W compartment due to pairings of Xm, and so on….
Thus:
Where sigma is the fecundity of the X females, and rho is the fecundity of the W females.
A model similar to this one can be found here (warning: there are many typos in their model formulation).
In the wolbach.R script, I run the simulation for 90 days (about the length of a rainy season), and see how many X’s and W’s remain in the population. With this model, we can examine the effect of releasing just male Wolbachia into the wild population, or a mix of male and female. Running 100 realisations with an initial wild population of 100, and a Wolbachia population of 100 equally mixed between male and female, produces the following plot:
Questions: what are the assumptions of the model? How can we make it more realistic?