Role of the potential landscape on the single-file diffusion through channels
Transport of colloid particles through narrow channels is ubiquitous in cell biology as well as becoming increasingly important for microfluidic applications or targeted drug delivery. Membrane channels in cells are useful models for artificial designs because of their high efficiency, selectivity and robustness to external fluctuations. Here we model the passive channels that let cargo simply diffuse through them, affected by a potential profile along the way. Passive transporters achieve high levels of efficiency and specificity from binding interactions with the cargo inside the channel. This however leads to a paradox: why should channels which are so narrow that they are blocked by their cargo evolve to have binding regions for their cargo if that will effectively block them? Using Brownian dynamics simulations, we show that different potentials, notably symmetric, increase the flux through narrow passive channels – and investigate how shape and depth of potentials influence the flux. We find that there exist optimal depths for certain potential shapes and that it is most efficient to apply a small force over an extended region of the channel. On the other hand, having several spatially discrete binding pockets will not alter the flux significantly. We also explore the role of many-particle effects arising from pairwise particle interactions with their neighbours and demonstrate that the relative changes in flux can be accounted for by the kinetics of the absorption reaction at the end of the channel.
pacs:83.10.Mj, 83.10.Rs, 87.16.dp
Transport of macromolecules through nano-sized pores and narrow protein channels is essential for cell function alberts2008molecular () while also becoming increasingly important in microfluidic applications shah2008 (); Pagliara2013 () or to understand drug delivery.Sugano2010 () Channels in cell membranes are remarkable for their high efficiency, selectivity and robustness with respect to fluctuations of their environment hille2001ion () and come in two flavours. Active transporters move their cargo by using cellular energy, e.g. from hydrolysing adenosine triphosphate or by harvesting concentration gradients of cell metabolites across the membrane. On the other hand, passive transporters are driven by the growth of entropy of the system as they translocate their specific cargo. Initially thought of as molecular sieves that select via the pore size to let the ‘right’ cargo simply diffuse through the channel, it is now well established that passive transporters achieve high levels of efficiency and specificity from binding interactions with the cargo inside the channel. A well-characterised example is the bacterial channel Maltoporin, where oligosaccharide transport is facilitated by an extended binding region.Schirmer1995 () Although many more examples of this phenomenon have since been discovered using a plethora of methods (e.g. ex-situ crystallographic studies,Kasianowicz2006 () indirect measurements of ionic currents Hilty2001 () and molecular dynamics simulations,Jensen2002 ()) the exact details of the mechanisms of passive transporters are still poorly understood.Kolomeisky2007 ()
Our work is motivated by a seeming paradox that arises when one considers the flux through a narrow channel, such as Maltoporin, which prevents particles from overtaking each other. Increasing the binding affinity between the channel interior and the cargo will prolong the time each particle spends inside the channel, hence reducing the flux and effectively blocking it. Why then would channels evolve to have binding regions for the molecules they have evolved to translocate? In this paper, we combine the results of Brownian dynamics simulations with theoretical arguments to show how developing binding regions inside a channel can indeed increase flux through narrow channels. We model the binding regions using a variety of potential-energy landscapes along the channel and investigate the dependence of the flux through the channel on the shape and depth of these potentials.
We first consider a single freely diffusing particle to tune our Brownian dynamics simulations in the setting where an exact analytical solution for the transport exists, applying various tests to the simulation procedure to ensure its proper reflection of the physical situation. We then investigate single-file diffusion through a channel to analyse the dependence of particle flux on the shape and depth of applied potentials. Finally, we demonstrate that we can account for the relative changes in flux by considering the diffusion-limited reaction kinetics of the absorption in a scheme Dorsaz2010 () based on the osmotic pressure along the channel alone.
To investigate the dynamics of colloid particles diffusing freely or in a potential, we need to solve the Langevin or corresponding Fokker-Planck equation. In a many-body problem like the one considered here, solving these equations directly is completely unfeasible and we therefore turned to Brownian dynamics (BD) simulations Chen2004 () using the LAMMPS package.Plimpton1995 () This allowed us to efficiently compute solutions of the equations from an ensemble of simulated stochastic processes with the same initial conditions.
Ii Free Brownian particle
The one-dimensional motion of a Brownian particle is described by the Langevin equation , where is the particle mass, its velocity and is the viscous drag coefficient. We further assume white noise with , where is the intensity of the stochastic force, satisfying the fluctuation-dissipation theorem. The general solution for the root mean square displacement of the particle is Uhlenbeck1930 ():
with the velocity relaxation time . In the overdamped (or diffusive) regime, with , we recover the Einstein result: , where one defines the diffusion constant . On the other hand, in the inertial regime , the displacement grows linearly with time: .
Brownian dynamics simulation
To verify that LAMMPS yields particle trajectories with the right statistical properties, we first simulated freely diffusing spherical Brownian particles of different sizes at different temperatures. The goal was to identify the crossover from the inertial to the diffusive regime at different temperatures and to verify the fluctuation-dissipation theorem.
Integration of the Langevin equation and the application of thermostat conditions was done via the fix_langevin routine.Schneider1978 () The free particles were simulated in a box with periodic boundary conditions and were assigned initial velocities drawn from a uniform distribution for the given temperature. The viscous drag coefficient was computed using Stoke’s law for a spherical particle at low Reynolds number: where is the particle radius and is the fluid viscosity.
Figure 1 shows the average root mean squared displacement of the particles computed as the average of 100 simulated trajectories per particle. Since there is no energy scale for the free particle, we used natural units and in this case set the diameter of the particles to and their mass density to that of water, yielding a mass of . We found that the statistical relative errors for the average trajectories were negligible for this number of simulations.
The crossover time was determined by inspection from the graphs. We can read off a crossing-over time of for particle #3 where the inertial response has fully died down; this is on the same order as the characteristic time scale .
We simulated particles of different sizes () at different temperatures (). Here it may be necessary to account for the fluid viscosity variation with temperature, and we used the empirical formula for water Al-Shemmeri2012 ():
We have confirmed that the particle trajectories generated by LAMMPS had the statistical properties expected from theory. We were also able to confirm that the crossover time is practically independent on the heat bath temperature, since the viscosity only depends weakly on temperature in the range that we covered in our simulations.
To verify the fluctuation-dissipation theorem, which is used to derive the Einstein relation , we computed the diffusion constant from a linear fit of the last three decades of each trajectory, i.e. for , obtained from 100 simulated particles using the SciPy Oliphant2007b () curve_fit routine. Figure 2 shows the product that was computed theoretically using the Stokes relation (see above) as a function of temperature, compared with the measured data. The predicted trend is observed, with a very small systematic offset that has been observed for a number of integration schemes in Brownian dissipative dynamics.Vattulainen2002 () Essentially, this is an artefact that arises from the coarse-graining of the microscopic properties of the fluid using a random force while imposing overall momentum conservation; this error is not significant for our purposes for two reasons: the algorithm produces trajectories in almost perfect agreement with the theory across a broad range of temperatures and furthermore, previous studies have shown that the effects due to integrator artefacts are only significant when the conservative forces of interest are comparable to the thermal fluctuations Vattulainen2002 () – all the potentials we apply will exceed energies of a few .
Iii Free particle in a channel
Having established the dynamics of Brownian particles and verified that LAMMPS generates trajectories with the desired statistics, we now turn our attention to the diffusion of free particles confined in a narrow channel.
We considered particles of diameter freely diffusing through a channel of radius . Note that from here on, we will use Lennard-Jones (LJ) units which render all quantities dimensionless by assuming that particle interactions follow the standard Lennard-Jones potential and setting the particle mass , the Boltzmann constant and and as defined above equal to 1 ljunits (). All masses are to be understood as multiples of these fundamental values, while other variables can be transformed to a dimensionless form by a scaling with an appropriate combination of the above, e.g. for time: . For a full list of conversion formulae see.ljunits () LJ units are widely used in computational physics and offer the advantage of being able to treat systems of different size and energy scales in one framework.
In our simulations, particles are modelled as spheres with a Lennard-Jones 12/6 type repulsion and no attractive interaction tail, that is, the LJ potential of pair interaction is truncated at the point of its minimum, . The channel radius is too small to allow particles to overtake each other, thus producing the single-file diffusion and reflecting the experimental fact that many metabolites will completely block their channels during the transport due to their tight fit.Nestorovich2002 ()
Figure 3(a) shows a 2D-projection of the setup of the simulation box (all simulations were carried out in full three dimensions). The ‘channel’ of length is the white area in the middle and is aligned along the -axis. Its walls interact with particles using the repulsive part of the Lennard-Jones potential, same as described above, thereby ‘softly’ preventing contact. Particles are inserted in the insertion region to the left (blue) if there is enough space. They then diffuse inside the cylinder. Once they have crossed the channel and entered the removal region to the right of the channel (blue), they are removed from the simulation. Underneath, in Fig. 3(b) is a plot of two example potentials and their corresponding force landscapes, to scale. To investigate the first passage time distribution, no potential was applied.
Distribution of first passage times
To check the physics of our channel setup, we looked at the distribution of first passage times, , of particles freely diffusing through the channel, i.e. without any applied potential. There is a classical result due to Lord Kelvin, obtained by the methods of images Lappala2013 () – however, it is only applicable when the particle is free to diffuse as far as necessary to the left of , while the first passage time is being tested by arriving at to the right of its entry point, see Fig. 3(a). In our case the passage is blocked to the left, so to find the probability for a particle one needs to solve the one-dimensional free diffusion equation with the boundary conditions: reflective wall, at , absorbing wall, at , and the initial condition for insertion: . The solution is:
where the constant is determined by normalisation. The survival probability for the particle to remain anywhere between 0 and , having started at is obtained by integration: . Given the boundary conditions, does not depend on anything happening outside the interval. Given the definition of the survival probability, the fraction of particles equal to is absorbed between and . This means that is actually the probability density of the time that takes the particle to reach for the first time. This distribution function is plotted in Fig.4, and it gives average first passage time .
We measured the passage times in 500 simulations of a single particle diffusion (to ensure no pair-interaction events could take place), where we inserted a particle in the insertion area and allowed it to freely diffuse to the end of the channel, the time for which we measured. We fit the distribution of first passage times to the resulting histogram of particle travel times in Fig.4, where only the diffusion constant and a normalisation act as fitting parameters. The data are in good agreement with the model, and we find a bare diffusion constant .
Interestingly, we found that when we considered a similar experiment where we inserted several particles into the channel one by one at a certain low frequency, the distribution of first passage times severely deviated from the free-diffusion result even at concentrations of just 2-3 particles inside a channel of length at any one time. This shows that many-particle effects caused by particle-particle interactions cannot realistically be ignored even at the lowest of concentrations, a point to which we will return at the end of this paper. In this particular case, the average first passage time was significantly increased at these concentrations, suggesting a smaller diffusion constant or higher effective resistance.
Iv Potential along the channel
Having established that our simulation setup produces physically meaningful results, we now turn to the dependence of the flux through a passive channel on the potential landscape inside it. We therefore made a series of experiments, in each of which we simulated the insertion of 100 particles in the channel at a rate . Since this is now a genuine multi-particle problem, simulation time increases accordingly. We therefore made use of the parallel computing capacities of LAMMPS.
Particles were only inserted if there was enough space in the insertion region; if a particle could not be inserted due to crowding at the channel entry, the insertion was skipped and the next attempt was made after a time interval of . Inside the channel, one out of a number of different potential shapes was applied with potential depth between in LJ units. The shapes of the potential are shown in the insets of Fig. 6. Note that potential ‘steps’ are modelled using functions, hence the names ‘Double tanh’ etc. Potential shapes include ‘continuous’ potentials, where the channel is modelled having a homogeneous attractive interaction along its entire length (single/double tanh, triangular potentials) or ‘discrete’ potentials, where the channel provides a number of discrete, spatially well-defined binding pockets, modelled as a Gaussian with a depth and a standard deviation of .
Despite the binding pockets being narrow, we can clearly see particle trapping occurring by looking at individual particle trajectories such as the one shown in Fig. 5, where displacement of a single particle along the channel is plotted on the x-axis in red, with time on the y-axis. The applied potential, a series of four spatially discrete binding pockets, each modelled as a Gaussians, is plotted in blue. We can clearly see that the particle is trapped by every binding pocket, spending most time in the second pocket from the left. However, the continuous insertion of particles to the left of the channel leads to a consistent movement of the tagged particle to the right.
We measured the cumulative number of particles that have crossed the channel as a function of time. An example of these measurements is shown in Fig. 6, where the cumulative number of translocated particles is plotted as a function of time for different potential shapes, all with the same . Insets show the exact form of the applied potentials for each experiment.
We can already make a number of observations from Fig. 6. First of all, a simple potential drop at the entrance of the channel leads to the expected acceleration and hence significantly increased the flux (red). Symmetrical potentials produce an increased flux compared to a channel with no potential () or different numbers of discrete binding pockets, which surprisingly produce a flux on par with the case. The number of translocated particles saturates at a value lower than 100, the total number of particles inserted during the simulation, because as particle insertion stops, exiting the channel becomes increasingly harder for the remaining particles since the pressure inside the channel decreases. The fact that the translocation number for a triangular potential saturates at a higher value than for the double-tanh potential supports this interpretation, since the double-tanh potential has a steeper wall at the end of the channel which will effectively block the channel.
Let us now discuss how the flux through the channel depends on the depth of the potential for the different potential functions discussed so far. We therefore define flux as the slope of a fit to the linear portion of the translocation plots in Fig. 6:
and repeat the analysis above for different potential depths. The average flux was computed from five experiments for given potential shape and depth, that simulated the translocation of 100 particles each. Simulating this number of particles stabilised the linear fits and resulted in the small relative errors that are plotted in Fig. 7 and enable us to refine some of the observations already made.
1. Symmetric potentials increase the flux of single-file diffusion. This is surprising at first sight, since the overall work done on the particle is zero. However, the symmetry of the potentials is broken by the pressure that the newly inserted particles exert on the particles near the end of the channel at the potential wall.Lappala2013 () It should be noted that this pressure emerges purely from the free diffusion inside the channel and has significant effects even at low colloid concentrations inside the channel. We are only inserting particles if there is free space at the beginning of the channel, as described above, so we are not actively pushing particles through the channel.
Furthermore, the flux through channels with symmetric potentials does not go below the flux through a free channel even for deep potentials.
2. The triangular potential profile outperforms the double-tanh potential. This is due to the fact that in the overdamped limit, after an impulse is exerted on a particle, it quickly relaxes back to normal diffusion. Effectively, Newton’s second law does not hold anymore and a small force over a longer time, pushing the dense region forward at the entry half of the channel, is more effective than a strong force over a short period of time.
3. The increase in flux with symmetric potentials is not due to some sort of Kramers-type barrier hoping. This is shown by the fact that the flux through a channel with a step at its end (‘tanh wall’) (see inset in Fig. 7) goes to zero for , where double-tanh and triangular potentials still outperform the case.
4. Narrow binding pockets do not alter the flux significantly, even though it is clear from individual particle trajectories that particles do get trapped in the binding pockets, as we see in Fig. 5.
Figure 8 gives the snapshots of final simulation frames to illustrate what is an ‘equilibrium’ situation in each potential profile . It shows that for a sufficiently deep attractive potential well, particles are retained in such a well, while the particles facing weaker binding forces escape and diffuse out of the channel. The final number of retained particles explains why the plateaus of different curves in Fig. 6 are below 100. These snapshots also help understand why the flux increases with the depth of continuous potentials (double-tanh or triangular). The process is analogous to the enzymatic action: although the energy barrier at the end of channel is prohibitively high (as illustrated by the complete vanishing of diffusive flux for the ‘tanh-wall’ potential in Fig. 7), when particles are confined at a high density in front of such a wall – they are forced to escape, pushed by the neighbours from the left.
It is also interesting to observe that at a constant temperature of our heat bath, when these potentials become excessively deep, the channel does get blocked: this occurs at for the double-tanh potential, and has to be inferred to occur at a much greater depth for the triangular potential, see Fig. 7.
V Theoretical considerations
We would like to gain a better insight into the data shown in Fig. 7 from a theoretical perspective. Since the basic features of single-file diffusion have already been discussed extensively elsewhere,Levitt1973 (); Kollmann2003 (); Burada2009 () we will focus our discussion here on the relative change in flux through a channel when we apply a potential.
It turns out that an effective way to pose this problem is to describe translocation as a reaction where the colloidal particles are absorbed by a “trap” , i.e. the channel exit, upon encounter. The problem of finding the flux through the channel becomes the problem of computing the rate of this reaction in a crowded single-file environment with applied potentials. For systems with spherical symmetry in the limit of infinitely diluted reactants , this is a classical problem of diffusion-controlled reaction kinetics, which was solved exactly by Smoluchowski,Smoluchowski1917 () producing the rate where is the density profile of reactants around the trap reaching the value at infinity.
In general, the reaction dynamics is governed by the diffusive Fokker-Planck equation:
which takes the many-body effects into account through the inhomogeneous density profile along the channel, which generates an osmotic pressure that acts as to spread the density profile via a force per particle .Zaccone2012 () These two parameters are non-linearly coupled via the collective diffusion coefficient , which makes an exact solution of this problem a formidable task even for numerics and makes us look for reasonable approximations.Zaccone2013 ()
In our 1-dimensional case, we take as the number of colloid particles per unit length along the channel at a given time (hence always for narrow channels), and we initially ignore the force due to the applied potential. With increasing , it takes a given particle longer to reach the channel exit, but once it is in the vicinity of the exit, its chance of reaching it and escaping increases due to the added gradient of osmotic pressure. Using this model for the analysis of their simulations, Dorsaz and co-workers Dorsaz2010 () showed that the reaction rate can be approximated well by the following expression:
where is the density at the beginning of the channel, and is the density a characteristic ‘encounter distance’ from the channel end at which the density of colloid particles acquires structure due to interactions (in other words, where the ideal-gas linear relationship stops being valid), see Fig. 9. Equation (6) has since been derived from first principles by Zaccone,Zaccone2013 () who finds a prefactor of instead of , but notes that that the two prefactors have the same dependence on which would indicate that the solution is qualitatively correct.
We write the osmotic pressure as a virial expansion in the density along the channel:
and compute from 1000 randomly selected snapshots of the simulations after the flux has equilibrated to its steady-state value. The virial coefficients and account for two- and three-body interactions between the particles in the channel which captures the essential dynamics since in the effectively 1-D system of the channel, the motion of a particle is dependent on the particle in front and the particle behind it.Lappala2013 () and were computed for a Lennard-Jones potential with that was used to model particle-particle interactions as described in section III. While ; computation of is more involved, but values are available.Bird1950 ()
A plot of the reaction rates computed from (6) is shown in Fig. 10. All values are normalised with respect to the reaction rate computed for no potential, . It is clear from the graph that these rates correctly predict the trends seen in the flux from the simulations (Fig. 7): there is no significant flux change with discrete pockets but a considerable increase with continuous potentials; the double tanh potential performs best at small potential depths while the triangular potential trumps at higher values of . The numerical range of the relative changes is good although it is systemically low by . This is a remarkable agreement given that at no point we explicitly introduced the form of the potentials and evaluate only at two discrete points, i.e. the beginning of the channel and very close to its exit. This shows that all the information about many-particle effects and the channel translocation with an applied potential is encoded in the steady-state density distribution, which in turn is controlled by the two virial coefficients and (or the two values and sampled near the beginning and near the end of the channel).
Our result that continuous, symmetrical potentials increase the flux significantly confirms earlier speculation that membrane channels in cells are most likely to provide a “molecular slide” Schirmer1995 () by organising discrete binding sites in succession, since having them isolated one after the other would provide little to no increase in flux as shown. Furthermore, our results can offer guidance for the design of artificial channels in microfluidic applications, where improving flux is often important and clogging can be a problem Lappala2013 ().
Any theoretical description of particle translocation has to account for both the applied potential and the crowding inside the channel. We have shown that it is possible to account for the relative changes in flux by considering the kinetics of the “absorption reaction” of particles exiting the channel, thus mapping the many-body problem to a two-body-interaction where crowding is modelled by the osmotic pressure inside the channel, without knowledge of the applied potential. However, this is more of an explanation a posteriori since it requires knowledge of the density profile along the channel. Further theoretical work will therefore have to focus on the development of methods to calculate these distributions not just for periodic boundary conditionsKollmann2003 (), but for more realistic geometries and boundary conditions in the presence of potentials in an attempt to predict particle flux without resorting to simulations.
We are grateful to Anna Lappala for many stimulating discussions. Simulations were funded by the Cavendish Laboratory teaching committee and performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England.
- (1) B. Alberts et al., Molecular Biology of the Cell, Garland Science, New York, 2008.
- (2) R. Shah et al., Mater. Today 11, 18 (2008).
- (3) S. Pagliara, C. Schwall, and U. F. Keyser, Adv. Mater. 25, 844 (2013).
- (4) K. Sugano et al., Nat. Rev. Drug Discov. 9, 597 (2010).
- (5) B. Hille, Ion channels of excitable membranes, Sinauer Sunderland, MA, 2001.
- (6) T. Schirmer, T. Keller, Y. Wang, and J. Rosenbusch, Science 267, 512 (1995).
- (7) J. J. Kasianowicz, T. L. Nguyen, and V. M. Stanford, Proc. Natl. Acad. Sci. 103, 11431 (2006).
- (8) C. Hilty and M. Winterhalter, Phys. Rev. Lett. 86, 5624 (2001).
- (9) M. O. Jensen, S. Park, E. Tajkhorshid, and K. Schulten, Proc. Natl. Acad. Sci. 99, 6731 (2002).
- (10) A. Kolomeisky, Phys. Rev. Lett. 98, 048105 (2007).
- (11) N. Dorsaz, C. De Michele, F. Piazza, P. De Los Rios, and G. Foffi, Phys. Rev. Lett. 105, 120601 (2010).
- (12) J. C. Chen and A. S. Kim, Adv. Colloid Interface Sci. 112, 159 (2004).
- (13) S. J. Plimpton, J. Comput. Phys. 117, 1 (1995).
- (14) G. Uhlenbeck and L. Ornstein, Phys. Rev. 36, 823 (1930).
- (15) T. Schneider and E. Stoll, Phys. Rev. B 17, 1302 (1978).
- (16) T. Al-Shemmeri, Engineering Fluid Mechnanics, Ventus Publishing ApS, 2012.
- (17) T. E. Oliphant, Comput. Sci. Eng. 9, 10 (2007).
- (18) I. Vattulainen, M. Karttunen, G. Besold, and J. M. Polson, J. Chem. Phys. 116, 3967 (2002).
- (19) The LAMMPS Team, Lammps documentation, http://lammps.sandia.gov/doc/units.html, 2014.
- (20) E. M. Nestorovich, C. Danelon, M. Winterhalter, and S. M. Bezrukov, Proc. Natl. Acad. Sci. 99, 9789 (2002).
- (21) A. Lappala, A. Zaccone, and E. M. Terentjev, Sci. Rep. 3, 3103 (2013).
- (22) D. Levitt, Phys. Rev. A 8, 3050 (1973).
- (23) M. Kollmann, Phys. Rev. Lett. 90, 180602 (2003).
- (24) P. S. Burada, P. Hänggi, F. Marchesoni, G. Schmid, and P. Talkner, Chemphyschem 10, 45 (2009).
- (25) M. Smoluchowski, Z. phys. Chem. 92, 129 (1917).
- (26) A. Zaccone and E. M. Terentjev, Phys. Rev. E 85, 061202 (2012).
- (27) A. Zaccone, J. Chem. Phys. 138, 186101 (2013).
- (28) R. B. Bird, E. L. Spotz, and J. O. Hirschfelder, J. Chem. Phys. 18, 1395 (1950).