Lecture 11: Revision of modelling approaches
Formulation of stochastic models
Remark: The gene regulatory network (toggle switch) demonstrated the importance of noise in biological systems. Without noise, no switch behaviour would be possible, and the system would always adopt only one state, starting from the same initial conditions.
Recap of Markov models
• Qi,j transition matrix: Entries are probabilities of a transition from state j to i
during time △t.
• The state space needs to be finite, which limits its application to the most simple
of models . Qi,j often take the form of αi,j △t or 1 − i̸=j αi,j △t.
• αk are called propensity functions where k denotes transitions from state j to i.
Example: Multiple independent and identically distributed (i.i.d.) channels instead of one
• Let Si = 0, . . . , N be i.i.d. state variables describing how many channels are open at a time
• The top left corner of a tri-banded Q matrix of size [N + 1] × [N + 1] would look like this (i, j refer to open states):
[1 − Nk+△t]{0 open → 0 open} [Nk+△t]{0 open → 1 open}
[k−△t]{1 open → 0 open}
[1 − (k− − (N − 1)k+)△]{1 open → 1 open} [(N − 1)k+△t]{1 open → 2 open}
.. .
• If N is large we immediately see △t needs to be small, because columns must add up to 1.
• A typical question is to find out how quickly a Markov chain converges to an equilibrium where π = Qπ, starting from a given intial distribution π0 of states.
• Noting that the solution will be the (right) eigenvector of Q associated with eigen- value λ = 1, the steady state can therefore be found by eigendecomposition of Q (see diffusion mapping), or simply by multiplying Q by itself until convergence is achieved.
• Remark: The left eigenvector associated with λ = 1 is always the 1 row vector consisting of ones only.
• The speed of convergence will be determined by the first eigenvector with |λk| < 1.
42
Example propensity functions
Examples for reactions, their associated propensity functions and changes in the state variables
k1
A→∅,α1 =A(t)k1 ,A(t+△t)=A(t)−1
k2
∅→A,α2 =k2,A(t+△t)=A(t)+1
k3
A+A→C,α3 =A(t)(A(t)−1)k3,A(t+△t)=A(t)−2,C(t+△t)=C(t)+1
k4
A+B → D, α4 = A(t)B(t)k4, A(t+△t) = A(t)−1 , B(t+△t) = B(t)−1,
D(t + △t) = D(t) + 1
(Remark: If we consider single molecules of A reacting with another molecule of A, we write A(t)(A(t) − 1),
which correctly results in 0 if there is only one molecule. In the law of mass action where concentrations are considered and we have large numbers of molecules we would simply write [A(t)]2).
Type
Model
Scheme
Prototype
Comments / Example
deterministic
ODE, chemical master equation
most generic formulation of reaction kinetics, can be time dependent
3 species with only up to 100 molecules: 1003 = 106 possible states and ODEs!
numerically approximations available, but only for selected problems
deterministic
discrete time Markov chain
use transition matrix Q, must meet Markov condition (no time dependence possible)
→− n →− P(t+n△t) = Q P(t)
For small problems where Markov chain converges to steady state distribution.
state space needs to be finite, single/multiple channel example
stochastic
discrete time, “naive” Monte Carlo
draw uniformly distributed random number r ∈ [0, 1], if
r < αi,j△t proceed with transition from state j to i.
△t must be very small, very inefficient
stochastic
exact, continuous time, Gillespie stochastic simulation algorithm (SSA)
draw two random numbers, one to predict time of next event, the other to pick next reaction
example of simple channel gating, and channel with diffusion (lab 4); different approximations of Gillespie’s method exist.
stochastic
discrete time, stochastic differential equation, approximation of noise through random forcing term, Langevin equation
Euler-Maruyama
dy(t) = f (y(t)) dt + σdW (t)
example of toggle-switch gene regulatory network
stochastic
random diffusion of particles,
Smoluchovski equation
X(t + △t) =
√
X(t) +
2D△tξx
examples of diffusion / subdiffusion (lab 4)
deterministic
discrete time, mean-field approach, replacing molecules with concentrations
integration using, Euler timestepping or higher order methods (Runge-Kutta)
dy(t)
= f(y(t) dt
Kai model for clock of cyanobacteria, Lotka-Volterra predator-prey model
Table 8: Modelling approaches
43
Lecture 12: Tracking problems in live cell microscopy
Examples:
• Population of dividing cells with one marker for cell nuclei, and additional markers for example for the cellular clock. Question: What is the relationship between cell division and the cellular clock?
• Kinetochore movements (centre part of chromosomes) during cell division. Ques- tion: What role play these movements in controlling errors in cell division?
• Molecular motors (kinesin) which attach to and walk along microtubules. Question: What can quantitative details reveal about the biochemical mechanisms powering this motion?
Virtual microscopy
• To extract time series data from moving objects we need efficient methods for object segmentation and feature extraction
• Simulated data is important to validate such methods
1. Simulate experimental data: To this end, hybrid models are often used, for example individual (agent) based models to simulate the cellular clock using a system of ODEs, the shape of cells or cell nuclei, together with movement using a stochastic models. To generate time series data, the noise and blur characteristics of a specific microscope are simulated.
2. Analysis of simulated/experimental data: The main steps comprise object detection, for example through template matching, and object tracking/feature extraction. Ideally, if the simulated experimental is representative of the true experimental data, the analysis will work on both equally well. The simulated data allows to increase real training data, for example when we are dealing with rare events like a few dividing cells amongst a majority of non-dividig cells.
Image formation in microscopy
• To simulate time series images, we need to understand some fundamentals of the image formation process first.
• The two main players are:
– Photon noise, which limits temporal resolution
– Blur introduced by diffraction at the aperture of the objective lens, which limits the spatial resolution
44
Photon noise
• In fluorescence microscopy, the release of photons from a fluorophore follows a Poisson statistics
• The probability distribution for counting p photons in an observation window of T seconds with ρ being the photon flux in photons per seconds and μ = ρT is:
P (p|μ) = (μ)p e−μ p!
• Define the signal to noise ratio (SNR) as
SNR = mean intensity of image or object
σ
with σ the standard deviation
• For Poisson noise the expected mean μ equals the variance (σ)2 so that
μ√ SNR=√μ= μ
given in dB: 10 log(SNR) dB
• Thus, SNR increases only slowly with increasing photon counts, that is more light, or longer integration time.
• Remark: In experiments with live cells, however, the amount of photons must be kept to a minimum as the introduced fluorescent markers and naturally occuring light absorbing cellular molecules can create toxic reactive oxygen species. High temporal resolution at low light levels will therefore result in significant amounts of noise.
Sampling Poisson noise
• To generate the appearance of a fluorescent object that follows a Poisson statistics we can sample from the distribution in the following way:
• For small μ, inverse transform sampling can be used P (p|μ) = (μ)p e−μ
p!
Algorithm: Poisson generator
(i)Set: p=0,q=e−μ,s=q
(ii) Draw uniformly distributed random number r ∈ [0, 1]
45
(iii) whiler>sdo{ p=p+1
q=qμ p
s = s + q %s: cumulative sum} return p
The Point Spread function (PSF)
• The PSF is an image of a point source, like a small fluorescent bead, in form of a blurred diffraction pattern.
• The PSF determines the ability to resolve two nearby points (maximum resolution: 200nm)
1 −x2+y2
• It is often approximated as a 2D or 3D Gaussian G(x, y) = 2πσ2 e 2σ2
• The blurring can be understood as filter operation, mathematically a convolution (denoted by the asterisk)
i=h∗o
where the object o is convolved with kernel h , the PSF, to yield the image i .
or in discrete form:
Creating a mock cell nucleus
i(x)= i(x, y) =
′ ′2′ h(x)o(x−x)d x
+m hp,q · ox−p,y−q p,q=−m
+m hp,q p,q=−m
• Create a ellipsoidal object, apply Poisson noise and convolution with a Gaussian (see example Matlab file)
• Intensity could be provided from a cellular clock SDE model for example
• Movement could be based on random diffusion (Smoluchovski)
• Overlap between nuclei could be prevented by employing a repulsive potential ∝ 1 r12
with r being the distance between the centre of nuclei (softshell with some possible overlap allowed)
% program to create a mock image of a cell nucleus
E = [40;60]; %axes of ellipse
shift = [100 ,100] ’; %shift centre
phi = pi/3; %rotation of ellipse
samples = 20; %discretisation of contour
I = 5; %intensity of nucleus
46
BG = 2; %intensity of background
theta = linspace(0,2∗pi ,samples); %angle parametrisation
R = E + 3∗randn(2,samples); %radius plus random fluctuation
x = [R(1,:).∗cos(theta);R(2,:).∗sin(theta)]; %generate parameterisation
x(:,end+1) = 0.5∗(x(:,1) + x(:,end)); %smooth endpoints M = [cos(phi) −sin(phi) %rotation matrix
sin ( phi ) cos ( phi ) ]
x =M∗x + shift; %apply rotation and shift
% generate a binary image and apply Poisson noise
% could use your own inverse transform sampling routine
out = uint16(BG + I∗poly2mask(x(1 ,:) ,x(2 ,:) ,256 ,256)); out = imnoise(out,’poisson ’); %create Poisson noise
% we use a 1D kernel for convolution and apply it in x and y
h=[0.129001 0.142521 0.151303 0.15435 0.151303 0.142521 0.129001] %sigma=5
out=conv2(out ,h’) %convolve in both directions out=conv2(out ,h)
imshow(out ,[0 max(out (:))]);
Simple tracking by template matching
• A simple but powerful method for template matching is based on correlation. A correlation is maximised if the template (kernel h) is a scaled version of the local intensity distribution in the image.
Convolution: i(x)conv = h ∗ o = +m hp · ox−p p=−m
object to be minimial.
+m +m
S = (hp −ox+p)2 = (h2p +o2x+p −2hp ·ox+p)
p=−m p=−m
+m +m +m
= (h2p)+ (o2x+p)− (2hp ·ox+p) p=−m p=−m p=−m
• The first term on the r.h.s. only depends on the size of the kernel, the second only on the object. The third is the correlation as by its definition. We see that the higher the correlation the smaller S will be.
• To detect point like objects like cell nuclei we can use a Gaussian kernel h which has a similar width as our object.
• Correlation with a discrete kernel works in pixel coordinates. To obtain a matching at subpixel resolution we can fit a continuous function like a Gaussian directly.
47
Correlation: i(x)corr = h ◦ o = +m p=−m
hp · ox+p
• If matching occurs, we can expect the squared distance between template and
• Typically, we would obtain the centres of the object through simple thresholding of the correlated image. This would serve as a starting point for a more refined segmentation, for example using graph based methods, clustering, active contour methods, etc.
• (Remark: Correlation could be wrongly suggested, simply because o is large. Also, for example, if the background takes high values, correlation will be high everywhere. Nor- malised cross-correlation (which is identical to computing Pearson correlation coefficient, not covered here) will overcome this problem.)
Tracking based on solving a linear assignment problem (LAP)
• The Hungarian or Kuhn-Munkres method is a polynomial time algorithm (O(n3)) to solve the optimal assignment problem
• The classical example is that of minimising the costs for assigning jobs to workers
• Asterisks denote the optimal assignment, with a minimum total cost of £490
Table 9: Assigning jobs to workers, minimising total cost
• We can use it to match objects in images between time frames τ and τ + 1;
• Multiple features can be integrated via a cost (or 1− similarity) function (distance between objects, intensity, size, shape, texture …)
Hungarian algorithm
Theorem 2. If a number is added to or subtracted from all of the entries of any one row or column of a cost matrix, then an optimal assignment for the resulting cost matrix is also an optimal assignment for the original cost matrix.
Algorithm follows https://www.mathworks.com/matlabcentral/fileexchange/20328-munkres- assignment-algorithm. A deeper understanding would require covering quite some back- ground on graph theory, for which, unfortunately, we don’t have time here.
• (0) Define matrix M of size n×n (n objects at time τ , n objects at time τ +1) with costs Mi,j for matching object i with j
• (i) Subtract the row minimum from each row
48
job 1
job 2
job 3
worker 1
£100
£120 *
£300
worker 2
£120
£130
£270 *
worker 3
£100 *
£130
£290
• (ii) Find a zero. If there are no starred (*) zeros in its column or row star the zero. Repeat for each zero.
• (iii) Cover each column with a starred zero. If all columns are covered then the matching is maximum.
• (iv) Find a noncovered zero and prime (’) it. If there is no starred zero in its row, goto Step v). Otherwise, cover this row and uncover the column containing the starred zero. Continue in this manner until there are no uncovered zeros left. Save the smallest uncovered value and goto Step vi).
• (v) Construct a series of alternating primed and starred zeros as follows: Let Z0 represent the uncovered primed zero found in Step iv). Let Z1 denote the starred zero in the column of Z0 (if any). Let Z2 denote the primed zero in the row of Z1 (there will always be one). Continue until the series terminates at a primed zero that has no starred zero in its column. Unstar each starred zero of the series, star each primed zero of the series, erase all primes and uncover every line in the matrix. Return to Step iii).
• vi) Add the minimum uncovered value to every element of each covered row, and subtract it from every element of each uncovered column. Return to Step 4 without altering any stars, primes, or covered lines.
See fully worked out example on module homepage
49