August 19, 2017, Saturday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Stochastic modeling of Pseudomonos syringae growth

From MathBio

Jump to: navigation, search


For this project, I attempted to recreate the results that were produced in the article entitled: "Stochastic Modeling Of Pseudomonos syringae growth in the phyllosphere"[1]

Biological Experiment

A mixture of Pseudoomonas syringae (often called P. syringae in paper) was sprayed onto 14 day old bean plants. The bacteria can be potentially harmful to plants as seen in the figure below. Part of the initial study was to better understand the growth of this bacteria on crop plants.This pathogen has been studied since the 1900’s and can ruin many plants such as: tomatoes, olives and bean plants. Part of the initial study was to better understand the growth of this bacteria on crop plants. This pathogen has been studied since the 1900's and can ruin many plants such as: tomatoes, olives and bean plants.


Mathematical Motivation

Figure 2

As the author quite plainly states, and with many mathematical modeling projects, "The aim of this work is to use a mathematical model to help understand(sic) the mechanisms" with in this case the mechanisms being growth of P. syringae on a leaf. Many other papers have been written on aggregates and cover a variety of different topics, from aggregates to social situations. One quite well known example would be “From Individuals to Aggregations: the Interplay between Behavior and Physics”[2] which studies and looks at a variety of aggregates on a both micro and macro scale. See Figure 2 for bluefin tuna school aggregate example. This paper not only chooses to investigate the aggregate formation but looks after examining basic results goes on to examine what happens as the initial distribution parameters change (due to nutrients) and as migration increases.

After looking around at other examples of bacterial growth, many authors may chose to use a probability distribution to model the data, but this article has aimed at understanding how it would arrive at such a distribution instead to understand the underlying mechanisms of the aggregate growth As the authors quite planely state their motive: "Rather than fitting a distribution to the experiemental data, we have opted for a dynamic approach that may allow us to elucidate the mechanisms that generate the observed behavior by Dulla and Lindow" While there may be some short comings in attempting to model this, it does allow for a much deeper understanding of how the aggregates came to be and many of the underlying biological actions that cause such a pathogen to resemble a specific distribution (in this case the [log normal distribution].

Mathematical Model

This model consists of three basic parts, the logistic birth-death process[3], the carrying capacity and the migration rate. A diagram of basic model is pictured below with no migration on top and migration on bottom.


Stochastic Model rather than Deterministic

The model they choose is stochastic, the authors explicitly state that they they will not discuss a deterministic model. Since all aggregates grow differently and independently, this creates a high degree of stochasticity which is not possible to do for many deterministic models. One precise reason they are using a stochastic process to model the experimental data is that you cannot say for certain when a cell will divide within a set amount of time, only that it has a certain probability of doing so.


  • Carrying capacity - the maximum colony size of each X_{i} colony. The carrying capacity was experiementally determined to be appoximately the log-normal distribution with parameters \mu =2.71 and \sigma =1.78 . Carrying capacity was randomly taken from Matlab with lognrnd(2.71,1.78) where as experiment used Randraw
  • Phyllosphere - the leaf surface or where the bacteria is inhabiting on the leaf
  • Aggregate - more than 1 cell colony
  • Monte Carlo Simulation - “computer algorithms that use random sampling to obtain numerical results”, was not used in my simulation but was in experimenters results[4]


  • Colony Migration (new colonies) do not affect old colonies
  • Cells can at any point migrate
  • Cells leave one by one
  • Colonies do no not effect one another, therefore each colony simulated independently
  • Two sources of new colonies are Inoculation (only initial) and Migration
  • N(t) at t=0, N(0) = 1
  • Sampling from lognormal distribution is unbiased. (Carrying Capacity has potential to be very very large even though experimentally this is impossible)


  • Growth/birth rate (.4): \alpha -
  • Death rate (.1): \delta
  • Migration: I
  • Lognormal parameter: \mu
  • Lognormal parameter: \sigma
  • Carrying Capacity: K
  • The number of colonies that were formed until time t: N(t) -
  • The number of cells in each colony: X_{i}(t)
  • Birth rate: \lambda _{i}, defined as \lambda _{i}=i\lambda (1-{\frac  {i}{K}})
  • Death rate: \delta _{i}, defined as \delta _{i}=i\delta


  • X(t) is a markov stochastic process with state spaces S = 0,1,...,K. i.e. memoryless with max population size K (integer step sizes)
  • Pr\left\{X(t+\Delta t)=j+1|X(t)=j\right\}=\lambda _{j}\Delta t - birth
  • Pr\left\{X(t+\Delta t)=j-1|X(t)=j\right\}=\delta _{j}\Delta t - death
  • Pr\left\{X(t+\Delta t)=k|X(t)=j\right\}=0 as \Delta t\rightarrow 0 Transition Probabilities Satisfy:
  • {\dot  p}_{x}^{{X_{i}}}(t)=p_{{x-1}}^{{X_{i}}}(t)\lambda _{{i,x-1}}+(\delta (x+1)+I)p_{{x+1}}^{{X_{i}}}-(\lambda _{x}+\delta x+I)p_{{x}}^{{X_{i}}}(t)
  • {\dot  p}_{m}^{{N}}(t)=p_{{m-1}}^{N}(t)(m-1)I-p_{m}^{N}(t)mI Solutions to this above equation can be obtained recursively and result in:
  • p_{n}^{{N}}(t)=e^{{-It}}(1-e^{{-It}})^{{n-1}}
  • Colony Formation Time: T_{i}=min\left\{t\leq 0;N(t)\leq i\right\}
  • X_{i}(t)=Y_{i}(t-T_{i}) if t>T_{i}, but 0 if t<T_{i}
  • Total Number of Colonies: N(t)=1+\sum \limits _{{i=1}}^{\infty }N_{i}(t-T_{i})
  • Y_{i} is birth-death process with Y_{i}(0)=1 initial population always 1
  • The birth or growth rate, is: \lambda _{{i,x}}, where \lambda _{{x,i}}=x\lambda (1-{\frac  {x}{K_{i}}}) since K_{i} is the carrying capacity for colony i, and x is colony size, as x \rightarrow K, growth or birth rate \lambda _{i} goes to 0. This allows for the colony size to asymptotically reach some carrying capacity K, as shown in figure 4
Figure 4

Simulation Results

Articles Simulation

  • The model was able to reflect the results achieved in the original experiment and is therefore useful for studying the mechanisms behind this pathogen and other similar ones.
  • The model reflects some interesting results such as, although the majority of colonies are in sizes less than 100, the majority of the populations is in colonies greater than 100.
    • This means that, there is a lot of population in only a few aggregates.
  • Integer values were used for carrying capacity

Simulation vs Experimental

  • In both experimental and simulation results:
    • Strong right-hand-skewed frequency distribution.
    • Large aggregates not frequent but account for majority of cells (more than 50% in colonies larger than 100)
    • Majority of colonies are small (colony size <50)
Simulation Results on Left and Experimental Results on Right

My Results

At first, I was able to recreate simulation in Simbiology, and with deterministic results, extremely similar to articles stochastic results except no variability (no birth and death fluctuation, population would reach carrying capacity and become flat). After finally being able to get stochastic results to work with a carrying capacity in Simbiology, I was able to very similarly recreate the results one of the authors showed me and my histograms were very similar to theirs as well with both deterministic and stochastic results.

My Results
My histogram on top and theirs on bottom with measurements taken at 4 different times

Matlab Code for my Simulation Results

% Create SimBiology Model.
m1 = sbiomodel('untitled');

% Add component(s) to SimBiology Model.
c1 = addcompartment(m1, 'unnamed');
s1 = addspecies(c1, 'N');
r1 = addreaction(m1, 'N -> 2 N');

% Add kinetic law object to the reaction.
k1 = addkineticlaw(r1, 'MassAction');

% Add component(s) to SimBiology Model.
p1 = addparameter(k1, 'k1', 1.0);

% Configure properties.
set(k1, 'ParameterVariableNames', {'k1'});

% Add component(s) to SimBiology Model.
r2 = addreaction(m1, '2 N -> null');

% Add kinetic law object to the reaction.
k2 = addkineticlaw(r2, 'MassAction');

% Add component(s) to SimBiology Model.
p2 = addparameter(k2, 'k2', 1.0);

% Configure properties.
set(k2, 'ParameterVariableNames', {'k2'});
set(s1, 'InitialAmount', 1.0);
set(p2, 'Value', 0.1);
set(p1, 'Value', .4);

% Initialize configset for analysis run.
cs = getconfigset(m1, 'default');

% RUNSIMULATION simulate SimBiology model, m1.

% Run simulation.
data = sbiosimulate(m1, cs, [], []);

% Assign variable name to object.
cs1 = getconfigset(m1, 'default');

% Configure properties.
set(cs1, 'SolverType', 'ssa');

% Initialize configset for analysis run.
cs = getconfigset(m1, 'default');

% Configure properties.
% time
set(cs1, 'StopTime', 50.0);
set(cs1, 'TimeUnits', 'hour');
% Run simulation.
data = sbiosimulate(m1, cs, [], []);

% Configure properties.
set(p1, 'Value', 10.0);

% Initialize configset for analysis run.
cs = getconfigset(m1, 'default');

% RUNSIMULATION simulate SimBiology model, m1.

% Run simulation.
data = sbiosimulate(m1, cs, [], []);
	%carrying capacity
% number of simulations:
num = 25;

for j=1:num
	mu = lognrnd(2.71,1.78);
	set(p2, 'Value', 0.1/mu);
	[t X] = sbiosimulate(m1);
	plot(t, X)
	hold on
axis([0 50 0 1800])

Suggestions for Model

While the model went on to further study Starvation driven migration, in which migration only occurred in colonies with greater than 50 cells to (so that only colonies with presumably fewer nutrients available would cause migration), an interesting model would be to compute it so that migration can happen at all times but becomes much more frequent close to carrying capacity.

For the most part, cell colonies either died within first 2-3 steps (because of stochasticity) or reached carrying capacity and stabilized.

Is this how real plant colonies occur or is there more fluctuation (taking advantage of susceptibility in plant)? Step size was in 1 hour increments, is there a more realistic option


  • Without a carrying capacity, the results are useless since cell colonies will grow exponentially.
  • Deterministic results quite similar to stochastic results but deterministic results arrive at carrying capacity slower.
  • Global (migration) and local (nutrient supply) effects affect population dynamics.
  • Population process that resembles experimental data was successful.
  • Possible use of model is for minimizing pathogens onset.

References and Further Reading


Along with the references, I found many links to good resources on a variety of topics covered (although I did not use all of these in my simulations):

Matlab code that helped me understand how to implement stochastic process without Simbiology: [1]

A slideshow introduction to Markov Chain and Sampling, very interesting: [2]

Monte Carlo Methods using Matlab, very in-depth and covers some parameter estimation : [3]

Logistic Growth Models: [4]

"A Short Course and Introduction to Dynamical Systems in Biomathematics" with great information and example Matlab files at end: [5]