Gillespie algorithm mtDNA population dynamics

gillespie.txt

				RATE OR		STOICHIOMETRY	
		REACTIONS	HAZARD		(W)	(M)	PROBABILITIES	
Synthesis	W -> 2W		ks*W		+1	 0	ks*W/CH
		M -> 2M		ks*M		 0	+1	ks*M/CH

Degradation	W -> NULL	kd*W		-1	 0	kd*W/CH
		M -> NULL	kd*M		 0	-1	kd*M/CH
						__________
Combined Hazard (CH):	(ks+kd)*(W+M)

W: Total number of wild-type mtDNA molecules in a cell		CN = W + M
M: Total number of mutant mtDNA molecules in a cell		Mutation Load (%) = 100 * M / CN
NULL: The void.  Nothingness.
CH: The combined hazard of all reactions occurring

GILLESPIE ALGORITHM
* Initialise cell at time t = 0 with W & M
* Calculate hazards for each of the four reactions
* Calculate the combined hazard
* Calculate the probabilities of each reaction occurring
* Randomly choose which reaction occurs by rolling 
	4-sided dice, weighted by probabilities
* Update W & M according to stoichiometry of selected reaction
* Calculate the time to next reaction event Dt ~ exp(1/CH)
* Update time t -> t + Dt