• Lab 4 (31.1.2020, to be completed in lab 5, which will have one additional smaller problem):
Stochastic modelling of Calcium ions entering the cytoplasm through a single membrane channel
Note: Care has been taken to choose reasonable values for all parameters and dimensions, in accordance with experimental findings. However, be pre- pared that some adjustments might be necessary if we find that computationally we will run into difficulties, because of, for example, numbers of molecules be- coming too large, or molecules moving too fast. The main rationale behind this practical is to get to know two fundamental Monte Carlo stochastic modelling techniques, namely the exact stochastic simulation algorithm (Gillespie), and Brownian motion. We will also explore subdiffusion, which is encountered when diffusion takes place in a confined space.
1) Channel gating
Simulate the stochastic opening and closing of a single ion channel and verify theaverageopenandclosedtime(⟨τ⟩= 1;⟨τ⟩= 1)aswellastheir
O k− C k+
theoretically expected distributions. Using the Gillespie algorithm employ the
following model:
k+
C (closed) O (open)
k−
Pick a uniformly distributed random number r ∈ [0, 1] and at time t compute τ = 1 ln(1/r) to determine the timepoints t+τ when the channel switches
C k+
from being closed to open, and τ = 1
C
ln(1/r) to determine timepoints t + τ O
O k−
when the channel switches from open to closed. Experiment with different values
of k+ and k−, but to get started you can use k+ = 0.1/ms and k− = 0.3/ms . Make sure that all your results contain the correct units.
2) Gillespie for modelling Calcium entry
The Gillespie algorithm presented in the lecture is a truncated version of the “proper” SSA algorithm. The reason for that is that we were only dealing with one reaction at a time, either a transition from C to O or the reverse. If we have i = 1,…,N multiple reactions that can occur , the Gillespie algorithm takes the following form.
a) Pick two uniformly distributed random numbers r1, r2 ∈ [0, 1].
b) Determine the timepoint of the next reaction event by computing τ = 1 ln1.
α0 r1
α0 = Ni=1 αi is the sum of all propensity functions αi, which tell how likely
it is that each reaction i is taking place.
c) Determine which reaction is taking place by finding j such that
1 j−1 1 j r2≥ αi and r2< αi
1
α0 i=1
α0 i=1
Update the numbers of reactants and products for the j’th reaction. Con- tinue with step (a)
We now include the translocation of ions through the channel into the cyto- plasm, with rate constant km = 5 L (units: L for litre, mM for millimolar,
mM ×ms
ms for milliseconds). The assumption is that ions pass through the channel as
long as it is open, and that there is a vast excess of ions in the extracellular space. That means we only need to keep track of the number of intracellular ions M and consider [Mextra] = 1.2mM/L, the concentration of ions in the extracellular space, constant.
Therefore, consider the following reactions with state variables C ∈ {0, 1}
(channel open or closed), O ∈ {0, 1} (channel closed or open) and M = 0, . . . , Mmax
the number of molecules that have passed through the channel. Initially there
are no ions present in the cytoplasm (M(t = 0) = 0). k+
i) Channel opens C −→ O ; C(t+τ) = C(t)−1 ; O(t+τ) = O(t)+1; M(t + τ) = M(t); propensity function α1 = k+;
k−
ii) Channel closed O −→ C ; C(t+τ) = C(t)+1 ; O(t+τ) = O(t)−1;
M(t + τ) = M(t); propensity function α2 = k−; km
iii) Ion passes through channel: Mextra −→ M ; C(t+τ) = C(t) ; O(t+τ) = O(t); M(t + τ) = M(t) + 1; propensity function α3 = km[Mextra]O;
Produce simulation runs (until M reaches Mmax = 1000, or a sensible value which makes plotting not too difficult!) and plots which demonstrate when the channel is open or closed and how the number of intracellular ions, M, increases. We use a comparatively low value for km, for real channels you could expect 106 ions entering per second, or more.
3) Simulating intracellular diffusion of Calcium ions
For a simulation run in problem (2) keep track of the times ti=1,...,500 when
the first 500 ions pass through the channel. We are now considering 2D diffusion
of these ions according to the Smoluchovski equation with D = 5 × 10−6 cm2 as sec
diffusion coefficient. The time steps △t in the Smoluchovski equation will be ti+1 − ti. We consider two scenarios, a) and b).
a) Diffusion in infinite halfspace, bounded by the membrane which contains our channel. Simulate the release of ions at the position of the channel (x = 0,y = 0) and allow diffusion in the halfspace defined by y ≥ 0. To prevent diffusion across the boundary y = 0, use reflection of ions at the boundary, that is Y (t + △t) = Y (t) + √2D△tξy. Plot trajectories of all ions.
b) Simulate release of ions into a square cell defined by the boundaries yb = 0, yt = 100μm , xl = −50μm and xr = 50μm, and subsequent diffusion. As in a) implement reflective boundary conditions (If the distances covered during timesteps △t turn out to be too big, that is using reflection cannot keep the ions in the box, you will need to take smaller steps △t. Alternatively (easier option), expand the size of the bounding box). Again, plot trajectories of all ions.
c) Plot the mean square displacement for the ions in a) and b) against time. Discuss your results using the formula msd(t) = 2nDtα.
2