Kinetic Diffusion
Fig. 1.Illustration of discretized model of ion diffusion (blue spheres within blueedged cells) near a mechanosensitive channel of small conductance (yellowedged cells) embedded in a biological membrane (red).
There are many ways of describing diffusion in a biological system. Free diffusion models based on the Green's function of the diffusion equation and random walk algorithms have been used to describe diffusion of biomolecules over the length scale of organelles and entire cells, on the order of micrometers, while highlydetailed Brownian dynamics and molecular dynamics have been used in simulations of diffusion over molecular length scales on the order of Ångstroms. Between the cellular and molecular length scales exists a middle ground where diffusion occurs over an intermediate length scale that is too large for atomistic models to be computed inexpensively and yet not large enough that the local potential and diffusivity fields can be neglected.
For example, the dynamics of ions translocating through an ion channel is affected by the diffusive approach of the ions in a large neighborhood hundreds of Ångstroms around the channel as well as the Ångstromscale variations in channel geometry and the potential and diffusivity fields near and within the channel. Few methods that efficiently simulate systems on these intermediate length scales and the corresponding typical diffusion time scales incorporate such detailed parameters. Driven by the challenge of simulating an ion channel on these scales, we combined the high resolution descriptive power of molecular dynamics with the efficiency of a finiteelement calculation in a particlebased method to describe systems on these intermediate length and time scales.
In order to incorporate the Ångstromscale details of the system, we first perform a longtime molecular dynamics simulation of the channel implanted in a membrane and immersed in a large waterbox containing the diffusing particles. The potential and diffusivity fields, U(r) and D(r) respectively, are extracted by averaging over time the particle concentration and mean squared displacement per unit time. The geometry of the system is implicit in the potential field, an example of which is shown in Figure 5. The potential and diffusivity fields are input parameters of the Smoluchowski equation, which the diffusion model described below is based on.
In our diffusion model, the system is first discretized on an irregular grid that assigns higher densities of cells to regions which require a detailed description, such as the neighborhood of an ion channel, and lower densities of cells to regions that are more uniform, such as empty space far away from the channel. With a good cell distribution, this discretization scheme makes the method more computationally affordable. Diffusing particles diffuse in the system by hopping across the grid in a series of time steps. For each hop, a particle chooses the target grid cell stochastically using a transition probability calculated from the Smoluchowski equation, given below:
The Smoluchowski equation describes the time evolution of probability distribution given an initial distribution, a potential field U(r) and a diffusivity field D(r). The equation is discretized by integrating over the cell currently occupied by the particle. The result of this integration is a set of probability flow rates R_{ij} into the current cell I_{i} from its nearest neighbors I_{j}, where j is an element of the set N_{i} of indices labelling the nearest neighbors of I_{i}. These probability flow rates depend on the potentials and diffusivities assigned by the preobtained fields to the respective cells, as well as the geometric properties of the cells, as described in Figure 2a. For each i, R_{ii}, which is the retention rate of particles within I_{i}, is the negative of the sum over j of R_{ji}, so as to conserve the particle number.
Setting then R_{ij} = 0 for j that are not in N_{i}, the kinetics of a particle moving from one cell to a nearest neighbor can be described by the master equation
Fig. 2. Schematic of the diffusion model. (a) Probability flow rates between I_{i} (pink) and nearest neighbors I_{j} (blue) specified in the master equation depend on the volume V_{i} of I_{i}, the interfacial area A_{ij} between I_{i} and I_{j}, and the distance between the cell centers w_{i} and w_{j}. (b) Transition probability T_{ji} from I_{i} to I_{j} for a particle initially in I_{i} after time Δt. The transition probability is obtained by integrating the master equation restricted to a neighborhood within a radius 2r_{Δt} of I_{i}.
The next step is to calculate, given that the particle is initially in I_{i}, the probability distribution over the grid after one time step Δt has elapsed. Such a calculation can become very expensive on large grids. Fortunately, it can be approximated heuristically that the particle is effectively confined within a certain neighborhood, of size dependent on the time step and characterized by length r_{Δt}, of its initial position. Hence, one can restrict the master equation (by including only the relevant matrix and vector elements) to the cells within this neighborhood. For each i, R_{ii} must be recalculated to preserve particle conservation. The restricted master equation is then integrated over time Δt, assuming that the particle is initially in I_{i}, to give T_{ji}, the probability of transition from I_{i} to I_{j}, where I_{j} is a cell within the neighborhood of I_{i} but not necessarily a nearest neighbor. Figure 2b illustrates a possible transition in one time step Δt within the restricted neighborhood. Finally, setting T_{ij} = 0 for all j not in the neighborhood of I_{i}, we can write for T_{ij} an analogy to the master equation in R_{ij}:
We tested our diffusion model by using it to describe simple systems, consisting of particles diffusing between two concentric spheres. Reflective boundary conditions were imposed on the outer sphere and absorptive boundary conditions on the inner sphere. The simulations were run with both a flat potential as well as a harmonic potential centered on the centers of the spheres. For different parameter values, the number of surviving particles, for which the analytic behaviors are known, were tracked over time and compared to the analytic behaviors. Figure 2 shows that for different values of time step size (2a) and grid densities (2b), the model agreed well with the analytic predictions for free diffusion. Figure 3 shows that the model holds up well after the addition of a harmonic potential, again for different values of time step size (3a) and grid densities (3b).
Fig. 3. Number of surviving particles for system with flat potential. Graphs for each case are superimposed and thus obscure each other.
Fig. 4. Number of surviving particles for system with harmonic potential. Graphs for each case are superimposed and thus obscure each other.
We also simulated the diffusion of K^{+} and glutamate (Glu^{}) ions through the mechanosensitive channel of small conductance (MscS) and measured conduction rates, with the same bias potentials used in BioMOCA simulations of the MscS in an earlier study. Our results are in agreement with the previous study in the case of K^{+}. In the case of the anion, since the previous study used Cl^{} as opposed to Glu^{}, there are significant differences in the conductance rates. Careful consideration of the reasons for these differences pointed to steric factors (Glu^{} is much bigger than Cl^{}) and more importantly, the lack of interparticle interactions in our model. Interparticle interactions will thus be the next step in the development of our model.
Fig. 5. (a) System used in molecular dynamics simulation. MscS (red) is embedded in a lipid bilayer (grey) and solvated in a waterbox. Also shown are potential maps of the same system for ions (b) K^{+} and (c) Glu^{}.
Publications

Publications Database A computational kinetic model of diffusion for molecular systems. Ivan Teo and Klaus Schulten. Journal of Chemical Physics, 139:121929, 2013. (15 pages).
Investigators
Related TCB Group Projects
Page created and maintained by Ivan Teo.