# Modelling Meso-Scale Diffusion Processes in Stochastic Fluid Bio-membranes

## Abstract

The space-time dynamics of rigid inhomogeneities (inclusions) free to move in a randomly fluctuating fluid bio-membrane is derived and numerically simulated as a function of the membrane shape changes. Both vertically placed (embedded) inclusions and horizontally placed (surface) inclusions are considered. The energetics of the membrane, as a two-dimensional (2D) meso-scale continuum sheet, is described by the Canham-Helfrich Hamiltonian, with the membrane height function treated as a stochastic process. The diffusion parameter of this process acts as the link coupling the membrane shape fluctuations to the kinematics of the inclusions. The latter is described via Ito stochastic differential equation. In addition to stochastic forces, the inclusions also experience membrane-induced deterministic forces. Our aim is to simulate the diffusion-driven aggregation of inclusions and show how the external inclusions arrive at the sites of the embedded inclusions. The model has potential use in such emerging fields as designing a targeted drug delivery system.   
PACS 87.20- Membrane biophysics.   
PACS 34.20- Interatomic and intermolecular potential.   
PACS 87.22BT- Membrane and subcellular physics.

Amphiphilic molecules, such as lipids and proteins, can self-assemble themselves into a variety of exotic structures in an aqueous environment . Computer-based simulation of the dynamics of these structures forms an interesting research area in the computational statistical mechanics of bio- material systems. One such structure is the phospholipid bilayer which represents the generic structure of all bio-membrane systems, both natural and artificial. These membranes can have thicknesses of only few nano-metres but linear sizes of up to tens of micro-metres, and can therefore be regarded as highly flexible, fluid-like, 2D continuum sheets embedded in a three-dimensional space.

Thermal fluctuations can induce shape fluctuations and shape transformations in membranes. For example, in the so-called budding transition  a spherical vesicle transforms into a prolate ellipsoid as the temperature increases. There is also the possibility that the spherical geometry becomes oblate, producing a shape similar to the biconcave rest shape of a red blood cell.

Bio-membranes regulate the recognition and transport processes to and from the interior of the cells, as well as between the cells, forming a barrier which all external particles arriving at the cell must cross. They contain a variety of integral (embedded) inhomogeneities, such as proteins and other macromolecules , that penetrate the thickness of the membrane and act as transport channels. We shall refer to these as the M-type inclusions. These inclusions are mobile and can freely diffuse across the membrane. Their presence, however, force the bilayer to adjust its thickness locally so as to match the thickness of the hydrophobic region of these inclusions , causing local deformations in the membrane geometry. The perturbations produced in the membrane shape due to these local deformations give rise to both short-range and long-range membrane-induced indirect forces between the inclusions. These forces act in addition to direct Van der Waals and electrostatic forces between the inclusions. The Long-range forces originate from the perturbations associated with the long wavelength shape fluctuations , whereas the short-range forces are associated with the local deformations in the immediate vicinity of the inclusions . These membrane-induced forces between the inclusions play a far more significant role, viz the direct molecular interactions, when the length scales involved are comparable to the size of the membrane. In addition to this mode of deformation, a membrane can also deform due to the tension at the amphiphilic molecules-water interfaces. This tension results in a change in the overall surface area of the membrane. A third mode of deformation also exists and this is associated with the bending elastic-curvature property of the membrane which distinguishes it from a sheet of simple fluid dominated by surface tension . Accordingly, two models to study the inclusion-induced local deformations have been developed. In the first model, the membrane energy is taken to consist of a contribution from the molecular expansion/compression due to the change in the thickness at the inclusion boundary, and also a contribution from the overall change in the surface area. Using this model, it is shown that  the inclusion- induced deformations cause exponential decays in the thickness of the membrane, extending from the inclusion-imposed value to the equilibrium thickness value, as shown schematically in Fig.1 for two rod-like inclusions. In the second model , the contribution of the membrane bending property is taken into account in the energy term, and it is found that this significantly affects the perturbation profile at the inclusion boundary as well as modifying the membrane- induced interactions. Evidently, an object supported by surface tension would have a different dynamics than one supported by the bending elasticity of the surface .

In addition to the M-type inclusions, membranes can also carry inclusions that lie on their surfaces  as shown schematically in Fig.2. These surface inclusions can represent objects that have arrived at the membrane from the outside and are therefore referred to as external inclusions. We shall refer to these as the S-type inclusions.

At the meso-scale, i.e. when the detailed molecular architecture of the membrane can be subsumed into a background 2D sheet, the free elastic energy of a symmetric membrane is described by the Canham-Helfrich Hamiltonian 

$$=d^2\sigma \sqrt{g}\{\sigma _0+2\kappa H^2+\overline{\kappa }K\},$$
(1)   
where

$`K=\text{det}(K_{ij})={\displaystyle \frac{1}{R_1R_2}},`$   
$`H={\displaystyle \frac{1}{2}}\text{Tr}(K_{ij})={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{R_1}}+{\displaystyle \frac{1}{R_2}}\right).`$ (2)   
are respectively the Gaussian and mean curvatures of the sheet, $`R_1`$ and $`R_2`$ are the two principle radii of curvature of the sheet, $`\sigma _0`$ is the surface tension, $`\kappa `$ is the bending rigidity, $`\overline{\kappa }`$ is the Gaussian rigidity, $`g`$ is determinant of the metric tensor and $`\sigma =(\sigma _1,\sigma _2)`$ is the 2D local coordinate on the sheet as opposed to the coordinates on the embedding space. The last term in (1) is, by Gauss-Bonett theorem, an invariant for closed surfaces implying that the dynamics of a membrane is not influenced by this term if its topology remains fixed. In what follows, we concentrate on membranes with fixed topology and drop this term. We then have

$$=d^2\sigma \sqrt{g}\{\sigma _0+2\kappa H^2\}.$$
(3)   
The study of a membrane whose free energy is described by (3) is facilitated by considering it to be nearly flat, i.e. its thickness to be much smaller than its linear size $`L`$. This is indeed what we mean by a meso-scale model of a membrane. We therefore take the membrane to be almost parallel to the $`(x_1,x_2)`$ plane, regarded as the reference plane. The position of a point on the membrane can then be described by a single-valued function $`h(x_1,x_2)`$ representing the height of that point. This simplification is achieved by writing the Hamiltonian (3) in the Monge representation  which gives for the mean curvature

$$2H=g^{3/2}[_1^2h(1+(_2h)^2)+_2^2h(1+(_1h)^2)2_1h_2h_1_2h],$$
(4)   
where $`_i\frac{}{x^i},i=1,2`$. We assume that the area of the membrane can fluctuate without constraint by setting $`\sigma _0=0`$ in (3). Consequently, using (4), the Hamiltonian (3) to leading order in derivatives of $`h`$ becomes

$$_0=\frac{\kappa }{2}d^2x(^2h)^2.$$
(5)   
This is the Canham-Helfrich Hamiltonian in Monge representation, expressed in terms of the height function of the membrane. It is the expression that we employ to describe the energetics of the membrane.

Employing a statistical mechanics based on (5) only, and ignoring the contributions from the expansion/compression and interfacial energies, the potential energy function

$$V_{MM}^T(R_{ij})=k_BT\frac{12A^2}{\pi ^2R_{ij}^4},$$
(6)   
was constructed  to describe the membrane-induced temperature- dependent long-range forces between a pair of disk shape M-type inclusions that can freely tilt with respect to each other. Another function was also constructed for long-range interaction between two S-type inclusions 

$$V_{SS}^T(R_{ij},\theta _i,\theta _j)=k_BT\frac{L_i^2L_j^2}{128R_{ij}^4}\mathrm{cos}^2[2(\theta _i+\theta _j)],$$
(7)   
where $`A=\pi r_0^2`$ is the area of an M-type inclusion of radius $`r_0`$, $`k_B`$ is the Boltzmann constant, $`R_{ij}`$ is distance between the centres of mass of two inclusions $`i`$ and $`j`$, $`L_i`$ and $`L_j`$ are the lengths of two S-type inclusions making the angles $`\theta _i`$ and $`\theta _j`$ respectively with the line joining their centres of mass (see Fig.2) and $`T`$ is the membrane temperature. It is evident that both of these membrane-induced potentials are attractive and fall off as $`R^4`$ with the distance. These expressions are derived for rod-like inclusions that are assumed to be much more rigid than the ambient membrane so that these inclusions can not move coherently with the membrane. The only degrees of freedom for the rods are rigid translations and rotations while they remain attached to the membrane.

So far, the modelling of bio-membrane dynamics decorated with inclusions has been mainly concerned with constructing potential energy functions such as those given in (6) and (7). An interesting problem, however, would be to use this information to simulate the space-time behaviour of inclusions in a membrane described by (5) and undergoing stochastic shape fluctuations. This is the problem that we address in this paper. This type of simulation can establish a direct link between the randomly changing membrane shape on the one hand and the inclusion dynamics on the other. In such a simulation, the thermodynamic phase behaviour of inclusions, such as their temperature-dependent aggregation, can be directly computed. This phase behaviour plays a crucial role in the functional specialisation of a membrane . Furthermore, information on the capture rate of the S-type inclusions, which could represent external drug particles, at the sites of the M-type inclusions can be obtained as a function of the changes in the environmental variables such as the ambient temperature. This type of meso-scale simulation when coupled with the Molecular Dynamics (MD) simulation of membrane patches near the inclusions at the nano-scale , can produce a seamless multi-scale model of the entire environment for many bio-molecular processes, starting with the arrival of external inclusions at the cell, their diffusion in the membrane, and finally their molecular docking at the site of the embedded inclusions.

To proceed, let us consider a 2D bio-membrane described by (5) containing both the M-type and the S-type inclusions. To make the membrane a stochastically fluctuating medium, we treat the height function in (5) as a stochastic Wiener process with a Gaussian distribution, whose mean and variance can be written as 

$$h(x_1,x_2;t)=0,$$
(8)   
$$h(x_1,x_2;t)h(x_1,x_2;t)=[h(x_1,x_2;t)h(x_1,x_2;t)]^2=2Dt$$
(9)   
where $`D`$ is the diffusion constant associated with the height fluctuations at the local position $`(x_1,x_2)`$ and represents the measure with which random fluctuations propagates in the local geometry. Such random height changes would cause a roughening of the membrane surface on molecular scales, and this has been observed in NMR experiments .

Assuming that this is the only stochastic process present in the membrane, it is then reasonable to suppose that this stochastic dynamics is communicated to the inclusions as well, and that their ensuing random motions are contingent only on these fluctuations. This implies that the mathematical point representing the centre of mass of an inclusion coinciding with the membrane point $`(x_1,x_2)`$ would also experience the same fluctuations and would diffuse with the same diffusion constant. To derive an expression for $`D`$, based on (5), we start with the static height-height correlation function obtained from (5). This is given by 

$$h(𝐪;0)h^{}(𝐪^{};0)=\frac{k_BT}{\kappa q^4}(2\pi )^2\delta (𝐪𝐪^{}),$$
(10)   
where $`\mathrm{}`$ is the thermal averaging with respect to the Boltzmann weight, $`\text{exp}(_0/k_BT)`$, and $`𝐪`$ is the wave vector of magnitude $`q`$. The corresponding dynamic correlation function can be obtained  by writing

$$h(𝐪;t)=h(𝐪;0)e^{\gamma _0(𝐪)t},$$
(11)   
giving

$$h(𝐪;t)h^{}(𝐪^{};t)=\frac{k_BT}{\kappa q^4}e^{2\gamma _0(q)t}(2\pi )^2\delta (𝐪𝐪^{}),$$
(12)   
where the damping rate, $`\gamma _0(q)`$, reflecting the long-range character of the hydrodynamic damping, is defined as

$$\gamma _0(q)=\kappa q^3/4\eta ,$$
(13)   
and $`\eta `$ denotes the coefficient of viscosity of the fluid membrane. In real space, (12) transforms  to

$$h(x_1,x_2;t)h(x_1,x_2;t)=\frac{k_BT}{4\pi \kappa }L^2e^{2\gamma _0t},$$
(14)   
where $`L`$ is the length of the membrane. This is the equal-time correlation function for membrane fluctuations. A similar model of an active fluctuating membrane in which the vertical displacements of the membrane satisfy a Langevin equation in the $`𝐪`$ space has also been proposed , and is its shown that a term similar to the static version of (12) contributes to the correlation function which also contains a contribution from non-equilibrium fluctuations. The latter is in the form of a $`q^5`$ term which dominates at long distances. Comparison of (9) and (14) yields the desired result

$$D=\frac{\left(\frac{k_BT}{4\pi \kappa }\right)L^2e^{2\gamma _0t}}{2t}.$$
(15)   
It should be emphasised that the association of a diffusive process with the membrane height function, and the resulting diffusion constant, is not analogous to the usual model of a diffusion process in which, for example, a particle diffuses through a medium, such as a fluid. Rather, what is suggested here is that the magnitude of a mathematical function representing the height of a mathematical point in the membrane is subject to random stochastic variations, and the diffusion constant is a measure of this variation.

When the M-type inclusions are present they produce exponentially decaying local deformations in the membrane geometry (see Fig.1). Correspondingly, the correlation function can be modified by a multiplicative exponential factor to

$$h(x_1,x_2;t)h(x_1,x_2;t)_{NI}=e^{r_0/R}h(x_1,x_2;t)h(x_1,x_2;t),$$
(16)   
where $`r_0`$ is the radius of an M-type inclusion, $`R+r_0`$ is the radius of the circular region around the inclusion with its centre coinciding with that of the inclusion, and $`NI`$ stands for near inclusion. It is evident that outside this region the exponential decay of the profile is negligible. Accordingly, within this circular region of radius $`R+r_0`$, the diffusion constant is also modified to

$$D_m=De^{r_0/R}.$$
(17)   
This equation implies that when the centre of mass an S-type inclusion enters a circular region of radius $`R+r_0`$ its diffusion coefficient goes over to $`D_m`$ and progressively approaches zero as the boundary of an M-type inclusion is approached. We can ascertain, as a first approximation, that this is how an M-type inclusion interacts with an S-type inclusion.

In our simulations, the equations motion of both the S-type and the M-type inclusions are represented by the differential equation of the the Ito stochastic calculus 

$$d𝐫(t)=𝐀[𝐫(t),t]dt+D^{1/2}d𝐖(t).$$
(18)   
This equation describes the stochastic trajectory, $`𝐫(t)`$, of the centres of mass of the inclusions in terms of a dynamical variable of the inclusions, $`𝐀[𝐫(t),t]`$, which is referred to as the drift velocity, and a term, $`d𝐖(t)`$, which is a given Gaussian Wiener process with the mean and variance given by

$`d𝐖(t)=0`$   
(19)   
$`d𝐖_i(t)d𝐖_j(t)=2\delta _{ij}dt.`$   
Equation (18) applies to each dimension of the motion. The Ito equation predicts the increment in position, i.e. $`d𝐫(t)=𝐫(t+dt)𝐫(t)`$, for a meso-scale time interval $`dt`$ as a combination of a deterministic drift part, represented by $`𝐀[𝐫(t),t]`$, and a stochastic diffusive part represented by $`D^{1/2}d𝐖(t)`$ and superimposed on this drift part. This equation resembles the ‘position’ Langevin equation describing the Brownian motion of a particle . The position Langevin equation corresponds to the long-time (diffusive time) configurational dynamics of a stochastic particle in which its momentum coordinates are in thermal equilibrium and hence have been removed from the equations of motion. Since we are interested in diffusive time scales as well, we can re-write (18) as

$$d𝐫(t)=\frac{D}{k_BT}𝐅(t)dt+D^{1/2}d𝐖(t),$$
(20)   
where $`𝐅(t)`$ is the instantaneous systematic force experienced by the $`i`$-th inclusion and is obtained from the inter-inclusion potentials, given in (6) and (7), according to

$$𝐅_i=\underset{j>i}{}_{𝐑_i}V(R_{ij}).$$
(21)   
We implemented (20) for our 2D simulations according to the iterative scheme

$`X(t+dt)=X(t)+{\displaystyle \frac{D}{k_BT}}F_X(t)dt+\sqrt{2Ddt}R_X^G`$   
(22)   
$`Y(t+dt)=Y(t)+{\displaystyle \frac{D}{k_BT}}F_Y(t)dt+\sqrt{2Ddt}R_Y^G`$   
where $`R_X^G`$ and $`R_Y^G`$ are standard random Gaussian variables chosen separately and independently for each inclusion according to the procedure given in , and $`F_X,F_Y`$ are the $`X`$ and $`Y`$ components of the force $`𝐅`$. For the S-type inclusions, we treated the angles in (7) as independent stochastic variables described by

$$\theta (t+dt)=\theta (t)+\frac{D}{k_BTL^2}\tau (t)dt+\frac{1}{L}\sqrt{2Ddt}\theta ^G,$$
(23)   
where $`\tau `$ is the torque experienced by an S-type inclusion and is given by

$$\tau _i=\underset{j>i}{}\frac{V^T(R_{ij},\theta _i,\theta _j)}{\theta _i},$$
(24)   
and $`\theta ^G`$ is the angular counterpart of $`R_X^G`$ and $`R_Y^G`$.

In the numerical simulations, recently reported in their broad outlines , we use a square membrane with $`L=40\mu `$m on its side. The other parameters used were set at $`\kappa =10^{19}`$ J and $`\eta =10^3`$ J sec m<sup>-3</sup> . These values correspond to the condition in which the bending mode of the membrane is important. From these data the damping coefficient, $`\gamma _0`$, in the real space of the membrane, can be obtained from (13). The simulation temperature was set at $`T=300^{}K`$, and the correlation (delay) time, $`t`$, over which the diffusion coefficient in (15) was calculated, was set at $`t=10^4`$sec. These data gave $`D=2.6\times \mathrm{\hspace{0.17em}10}^9`$ m<sup>2</sup>sec<sup>-1</sup>, which is in close agreement with the value of $`D4.4\times \mathrm{\hspace{0.17em}10}^9`$m<sup>2</sup>sec<sup>-1</sup> obtained at the molecular level via an MD simulation of a fully hydrated phospholipid dipalmitoylphosphatidylcholine (DPPC) bilayer diffusing in the $`z`$-direction . To justify our choice of the correlation time, $`t=10^4`$ sec, we recall that the time scale of a stochastic particle, $`t_D`$, of mass $`m`$ is usually determined from the relation 

$`t_D={\displaystyle \frac{mD}{k_BT}}`$ (25)   
Since $`t_D`$ is normally of the order of $`10^9`$sec, then for the criterion of long-time dynamics, employed in our model (cf (20)), to be justified the correlation (diffusive) time scale, $`t`$, in (15) has to satisfy the condition 

$`tt_D.`$ (26)   
For our calculated value of $`D`$ and our choice of the inclusion mass $`m=10^{12}\mu `$g, corresponding to an inclusion of length $`L_i=0.1\mu `$m, we obtained a value of $`t_D=0.6\times 10^9`$sec, showing that our choice of the correlation time was appropriate to satisfy the condition in (26). The radius of an M-type inclusion was set at $`r_0=0.01\mu m`$, and the inclusions were all equal in length.

The stochastic trajectories of the inclusions were obtained in a set of five simulations. The simulation time step, $`dt`$, in (22) was set at $`dt=10^9`$sec, and each simulation was performed for $`4\times \mathrm{\hspace{0.17em}10}^6`$ time steps, i.e for a mesoscopic interval of $`4000\mu `$sec. The total number of inclusions considered was $`36`$, consisting of $`13`$ S-type and $`23`$ M-type.

In the first simulation, we computed the random motions of the S-type inclusions in a membrane devoid of the M-type inclusions. This was done in order to observe the details of the drift-diffusion motion over mesoscopic time scales. Figure 3 shows the stochastic X-Y trajectories of a sample of 4 S-type inclusions plotted on a micron scale up to the end of the simulation time. In addition to the drift motions, represented by the second terms in (22), the random Brownian-type variations, emanating from the membrane shape fluctuations, are superimposed on this drift motion and are clearly visible over the mesoscopic length and time scales. Figure 4 shows the snapshots of a small patch of the membrane with both the S-type (white spheres) and the M-type (black spheres) inclusions. In this, and subsequent figures, the solid spheres refer to the centres of mass of the rod-like inclusions. In the initial state the outer M-type inclusions were regularly positioned, whereas the inner ones were randomly distributed. The S-type inclusions were all distributed completely at random. Figures 4a to 4c refer to the simulation in which the M-type inclusions were pinned to the membrane, i.e. were static, and only the S-type inclusions were mobile, and figures 4d to 4f refer to the simulation in which both the M-type and the S-type inclusions were mobile. The initial states in both simulations, figures 4a and 4d, were the same. The snapshots were obtained from dynamic simulations, akin to a MD simulation, covering the entire simulation time interval. These snapshots were recorded at $`2\times 10^3`$ sec intervals, with figures 4c and 4f referring to the final states reached at the conclusion of the simulations after $`4\times \mathrm{\hspace{0.17em}10}^6`$ time steps. The animation of a complete run showed clearly the stochastic motions of the inclusions, and how the S-type inclusions approached the M-type inclusions and were captured at the site of the M-type inclusions. An examination of figure 4 shows that for the case of dynamic M-type inclusions, a larger number of the S-type inclusions were captured at the M-type inclusion sites, i.e. the number was some 4 times higher than in the static case at the same temperature. We adopted the method of counting an S-type inclusion as a captured inclusion when its centre of mass coincided with that of an M-type inclusion. The numerical algorithm then transformed the colour code of that S-type inclusion from white to black. Figures 4e and 4d also show some diffusion-driven local aggregation of the M-type inclusions. The capture of the S-type inclusions can be viewed as the first stage in the molecular docking process which will eventually transfer these inclusions into the interior of the cell.

To examine the membrane response to temperature changes, two of the simulations were performed at different temperatures. Figure 5 shows the snapshots of these simulations at $`T=100^{}`$K (a to c) and at $`T=350^{}`$K (d to f). Figures 5d to 5f clearly show that both the number of captured inclusions and the aggregation of the M-type inclusions were affected by these temperature differences, as can be seen by comparing figures 5c and 5f.

To sum up, although many dynamical aspects of membrane-like surfaces have been addressed in the past , it is only relatively recently that attention has focused on the dynamics of membranes with inclusions. To our knowledge no computer-based simulation of this dynamics has been reported so far. In this paper we constructed a meso-scale model of a generic bio-membrane based on the Canham-Helfrich curvature-energy formalism. We treated the height function of the membrane as a stochastic Wiener process whose correlation function provides the relevant diffusion constant describing the membrane fluctuations. Two types of inclusions, one mimicking the internal embedded type and the other the external floating type, are carried by this membrane. These inclusions experience the same stochastic fluctuations as those experienced by the membrane itself, resulting in the transformation of their deterministic (drift) space-time dynamics into a stochastic Langevin-type dynamics described by the Ito stochastic calculus. A set of dynamic simulations, resembling the standard MD simulations, were performed to investigate the phase behaviour of these inclusions. In addition to stochastic forces, these inclusions also experience deterministic interactions described by inter-inclusion potentials varying as $`1/R^4`$ with their separations. The simulation results clearly indicate that the capture and aggregation rates change with the temperature and that the embedded mobile inclusions capture a greater number of the floating inclusions. A further extension of the present work would be to include the influence of the surface tension, as well as the bending rigidity, by keeping the corresponding term in (3). This will constrain the fluctuations in the surface area of the membrane and would have a direct bearing on the inclusion dynamics.   
The second author(HRS) is grateful to the UK Royal Society for financial support through a visiting research fellowship and to the School of Computing and Mathematical Sciences (Greenwich University) for their hospitality. Both authors acknowledge useful discussions with Professor E. Mansfield FRS on the dynamics of objects supported by surface tension.

Figure captions   
Figure 1: Two rod-like embedded (M-type) inclusions vertically placed in an amphiphilic fluid membrane. The inclusions impose exponentially decaying thickness-matching constraints on the bilayer at the inclusion boundary. Heavy solid lines represent amphiphilic molecules. Figure based on .   
Figure 2: Two rod-like surface (S-type) inclusions lying on the surface of the membrane. The rods have lengths $`L_1`$ and $`L_2`$, widths $`ϵ_1`$ and $`ϵ_2`$ and making angles $`\theta _1`$ and $`\theta _2`$ with the line joining their centres of mass. Figure based on .   
Figure 3: A small patch of the membrane showing the stochastic X-Y trajectories obtained from equation (22) for a sample of four S-type inclusions without the presence of the M-type inclusions. Both the drift and diffusion motions can be clearly distinguished.   
Figure 4: A set of snapshots, obtained from dynamic simulations, showing the capture of rod-like S-type inclusions (white spheres) at the rod-like M-type inclusion sites (black spheres) for static (a to c) and dynamic (d to f) M-type inclusions. The aggregation of the M-type inclusions can also be observed (d to f). Only the centres of mass of the inclusions are shown.   
Figure 5: A set of snapshots, obtained from dynamic simulations, showing the capture of rod-like S-type inclusions at the rod-like M-type inclusion sites for mobile M-type inclusions at $`T=100^{}`$K (a to c) and $`T=350^{}`$K (d to f). Only the centres of mass of the inclusions are shown.