Journal of Undergraduate Research
Volume 7, Issue 4 - March/April 2006

Molecular Modeling of Nanoparticle Transport across Lipid Bilayers

Ryan A. Tasseff and Dmitry I. Kopelevich

ABSTRACT

The rapid development of nanotechnology and increased ability to manufacture materials and devices with nanometer-feature size have led to exciting innovations in many areas, including within the medical and electronic fields. However, very little is known about the possible health and environmental impacts of manufactured nanomaterials. Recent experimental data suggest that some of the manufactured nanomaterials, such as fullerenes and their derivatives, are highly toxic even in small concentrations. The goal of the current work is to develop a theoretical and computational framework for fundamental understanding of the mechanisms responsible for the toxicity of nanomaterials. Here, we present the result of the initial stage of this investigation aimed at understanding the interactions between manufactured nanoparticles and cellular membranes. In the current work, we applied course-grained molecular dynamics simulations to investigate the mechanism of transport of carbon-based nanoparticles across a phospholipid bilayer. This transport rate is shown to significantly decrease with particle size, and the residence time of a nanoparticle of radius 0.88 nm is shown to be as large as 0.01 sec. This long residence of the particle within a bilayer may lead to the disruption of the membrane integrity as suggested by recent experimental evidence.

INTRODUCTION

Understanding the interactions between manufactured nanoparticles (NP) and biological cells is of increasing importance and may lead to many exciting innovations. For example, recent studies have shown that certain NP, including fullerenes, can serve as an effective means of drug delivery.1

Because of the growing number of applications, the production of nanomaterials is expected to significantly increase. But what dangers, if any, does mass production pose to the environment and general health? The truth is that the possible negative impact of these particles is not yet well understood. However, recent studies have revealed some alarming trends. Specifically, it has been observed that exposure to an aqueous solution of low concentration (0.5 ppm) of uncoated fullerenes leads to significant oxidative brain damage of large mouth bass within 48 hours.2 These results have increased the general concern over nanomaterials and have motivated several groups to work towards a more comprehensive understanding of their health and environmental impacts.

Further studies have shown that the presence of fullerenes, at concentrations as low as 20 ppb, induce membrane rupture and eventual death of human dermal fibroblast cells.3 The proposed cause of membrane rupture is the peroxidation reaction between reactive oxygen species and membrane lipid tail molecules. Although the reaction is not directly caused by fullerenes, it was shown that these NP do increase the concentration of reactive oxygen species in water. The apparent cytotoxicity of fullerenes is a sufficient motivating factor for better understanding the interactions between NP and living cells.

This study represents an initial step toward development of a fundamental understanding of the interaction of manufactured nanoparticles with cellular membranes. Here, we have investigated the transport mechanisms of NP through cell membranes. More specifically, we have focused on the transport of spherical, carbon-based NP across a model membrane, phospholipid bilayer. Transport of biomolecules into the cell interior usually occurs through a complex mechanism involving proteins. However, some molecules, such as water, are known to cross the membrane through direct, or basal, permeation.4 In the current study, we investigate the basal permeation of the NP through the lipid bilayer, as illustrated in Figure 1.

Figure 1. Schematic representation of a nanoparticle moving through a lipid bilayer

Figure 1. Schematic representation of a nanoparticle moving through a lipid bilayer. The z axis is perpendicular to the bilayer surface.

The goal of this work is to investigate effects of particle size on the transport rates and mechanisms of transport through the bilayer. The small length scales of interest in this work motivate the use of molecular dynamic (MD) simulations, which allow investigations of transport processes at a molecular level. It is anticipated that the time scales of NP transport are beyond the reach of direct MD simulation, and, therefore, it is necessary to employ an indirect method of calculating transport properties. To do this, we employed the constrained simulation techniques that were used earlier5,6 to investigate transport properties of smaller molecules, such as water, in similar bilayer systems.

These earlier studies employed a fully atomistic MD model to study the transport of small molecules such as water and oxygen. The results suggest that the transport of small molecules through the lipid bilayer is an activated process that takes place by hopping through the voids between the lipid molecules. It is unlikely, however, that larger particles will find voids of sufficient size for transport.7 Therefore, we anticipate a qualitatively different transport mechanism for the larger particles. Our hope is to obtain a more detailed understanding of the effects of particle size on the transport processes.

METHODS

Model Details

As discussed above, we investigate the nanoparticle transport using the molecular dynamics (MD) simulations. Our studies are performed using a coarse-grained molecular dynamics (CGMD) model, which approximates small groups of atoms as a single, united atom. Several such models have been introduced and applied to simulations of various complex molecular systems.8-12 The development of coarse-grained molecular models and the improvement of the quantitative agreement with full atomistic models and experimental data are subjects of active ongoing research.11, 13 In this work, we will use the model for water and phospholipids proposed by Marrink and Mark.12 Simulations with this model are shown to yield good agreement with experiments and atomisticlly detailed simulations for such quantities as density and elasticity of lipid bilayers.13

This model, on average, employs a four-to-one mapping to represent groups of atoms as a single reaction site, or bead. The following four main types of beads are considered: polar (P), nonpolar (N), apolar (C), and carged (Q). For example, four water molecules are represented as a single polar bead and four methyl groups of the lipid tail are represented by a single hydrophobic apolar bead. The interactions between non-bonded beads are modeled by the Lennard-Jones potentials (U),

Equation (1)

where the effective bead diameter (σ) is 0.47 nm, unless otherwise stated. The strength of interaction (ε) is limited to five possible values ranging from weak (1.8 kJ/mol) to strong (5 kJ/mol). Figure 1 shows a plot of Lennard-Jones potentials as a function of the distance, r, between the two interacting coarse-grained beads.


Figure 2. Interaction potential between two coarse-grained beads as a function of distance, r

Figure 2. Interaction potential between two coarse-grained beads as a function of distance, r. The solid line shows strong interactions between two hydrophilic beads (note the deeper well). The dashed (dotted) line shows intermediate (weak) interactions. In these plots _ = 0.47 nm

The large value of ε corresponds to strong attractions between hydrophilic beads, whereas the small values of ε are used to model interactions between the hydrophobic and hydrophilic beads that are dominated by repulsion. Finally, the interactions between two hydrophobic particles are modeled using an intermediate interaction strength (ε = 3.4 kJ/mol).

In addition, electrostatic interactions between charged beads are taken into account by the screened electrostatic potential,

Equation (2)

Here f = 138.9 kJ mol-1 nm e-2 is the electric conversion factor, qi is the screened charge of particle i, and εr is the relative dielectric constant.

Interactions between chemically bonded beads are modeled by harmonic potentials for the bond length and bond angle:

Equation(3)


Equation(4)

Where the force constant of eq. (3) is Kbond = 1250 kJ/(mol nm), r is the instantaneous bond length, rbond is the equilibrium bond length (usually 0.47 nm), the force constant of eq. (4) is Kangle = 25 kJ/(mol rad2), θ is the instantaneous bond angle and θ0 is the equilibrium bond angle.

Species Details

The cell membrane is modeled by a dipalmitoylphosphatidylcholine (DPPC) lipid bilayer. The DPPC molecule is comprised of the hydrophilic head group, which includes ester backbone, and two hydrophobic tails. The coarse-grained model for this molecule, shown in Figure 3a, contains two charged beads (Q) for the zwitterionic PC head group, two non-polar beads (N) for the glycerol ester backbone and four apolar beads (C) for each of the two tails. The lipid bilayer is surrounded on both sides by water. The coarse-grained model uses a single polar bead to represent four water molecules.

Coarse-grained model of the DPPC lipid molecule

Figure 3. (a) Coarse-grained model of the DPPC lipid molecule; (b) Initial random dispersion of DPPC beads in water; (c) Self-assembled DPPC lipid bilayer. For clarity, the water molecules are not shown.

A computational model for lipid bilayer is prepared through self-assembly. We performed MD simulations of a DPPC lipid solution in water in a simulation cell of size 10 x 10 x 10 nm3. The initial conditions for this simulation were chosen to be a random dispersion of lipids in water (see Fig. 3b). We found that lipids self-assemble into a bilayer shown in Fig. 3c within a few nanoseconds.

To facilitate the discussion to follow, we introduce the system of coordinates so that the z axis is perpendicular to the bilayer surface (see Fig. 1), and the x-y plane corresponds to the bilayer center of mass. In this report we present our results for spherical nanoparticles of various sizes. The goal of this work was to model carbon-based nanoparticles, such as fullerenes, and, therefore, we chose the model parameters, such as the range of effective diameters, to mimic the fullerenes.14 Therefore, nanoparticles were modeled as hydrophobic by assigning apolar interaction strengths, the same interaction strengths used for methyl groups. The nanoparticle effective diameter, σ, was varied between 0.47 nm and 1.5 nm.

Simulation Details

The CGMD simulations were performed using the Gromacs V3.1.1 software package.1515 The system contained 266 DPPC molecules in a bilayer configuration and 4007 water beads. The periodic boundary conditions in all dimensions were assumed. Temperature and pressure were maintained at 323 K and 1 bar using the Berendsen temperature and anisotropic pressure coupling methods.16 Anisotropic pressure coupling allowed maintaining zero surface tension of the bilayer membrane. The time step for the numerical integration was 0.04 ps.

Constraint Method

Transport of even small molecules across a lipid bilayer takes place on a time scale that is out of reach of MD simulations. It is therefore not practical to directly study the transport effects, and several techniques have been developed to overcome this obstacle. These techniques include umbrella sampling,17 coarse kinetic simulations,18,19 and the constrained mean force (CMF) method. In our investigation, we employed the CMF method. This method has been successfully applied to similar systems.5,6 The main idea of the method is to apply a constraint to force the NP of interest to remain at a given position, z0, along the reaction coordinate, z. The reaction coordinate is along the axis perpendicular to the bilayer surface and is shifted such that z = 0 corresponds to the bilayer center of mass (see Fig. 1).

From the constrained simulations, we obtained the force, F(z0, t), required to constrain the NP. Statistical analysis of F(z0, t) outlined in Section 2.5 enabled calculated of the parameters for a model describing the nanoparticle transport. The constraint is implemented using the SHAKE algorithm,20 a built-in function of the Gromacs simulation package.

Calculation Details

The main assumption underlying our analysis was that NP transport can be described by the generalized Lagenvin equation:21

(5)

Here, m is the mass of the nanoparticle, γ is the time-dependent friction coefficient, G is the free energy, and Γ is the normally distributed random force with zero mean and the autocorrelation function obeying the fluctuation-dissipation theorem,22

Equation (6)
Equation(7)

The brackets denote the ensemble average, kβ is the Boltzmann constant, and T is the system temperature. The fluctuation-dissipation theorem, eq. (7), provides a relationship between Γ and γ. As discussed in Friction Coefficient and Diffusivity, this relationship allows computation of the friction coefficient from the random force autocorrelation function, which is readily measurable from MD simulations.

Typically, the correlation time for Γ is very small (as discussed in section 3.2), and γ can be replaced by the time-independent friction coefficient, γ. This allows approximation of eq. (5) and eq. (7), respectively, as:

Equation (8)

Equation(9)

where δ is the Dirac delta function.

The computation of the terms in eq. (8) from the constraint MD simulations is discussed below.

Free Energy

If z is held constant by the application of a constraint force, F, then the first and second terms of eq. (8) are zero. Further, the ensemble average of the random force term is zero [see eq. (6)]. Therefore, the ensemble average of eq. (8) yields the following relationship between the constraint force <F> and the free energy:

Equation (10)

In analysis of our simulations, we replace the ensemble averaging by the time averaging using the ergodic assumption.

Friction Coefficient and Diffusivity

To fully reconstruct eq. (8), it is necessary to obtain the friction coefficient γ. To this end, we apply the relationship in eq. (7) to calculate the time-dependant friction coefficient γ from the autocorrelation function of the random force Γ,

Equation (11)

The time-independent friction coefficient γ is obtained by integrating eq. (11):

Equation(12)

Further, the obtained friction coefficient γ is directly related to diffusivity,21

Equation (13)

The random force Γ(z,t) can be obtained by computing the deviation of the constrained force from its mean value,

Equation (14)

Rate of Transport

Ultimately, we wish to consider the mean time required for the NP to move across a bilayer, τtrn. It can be shown21 that eq. (8) leads to the following expression for the average transport time τtrn between points z1 and z2:

Equation (15)

where z1 is the point of entrance, and z2 is some point of exit. In what follows, the integrals in this equation are calculated numerically. Eq. (15) is derived assuming constant diffusivity, D. The validity of this assumption is discussed in section 3.2.

Note that the expression, eq. (15), can be reduced to the more familiar Arrhenius relationship derived earlier by Kramers:23

Equation (16)

This latter result relies on several assumptions including: large separation of energies at the minimum and maximum of the free energy profile and parabolic shape of the free energy profile near the minimum and maximum. In our system the first statement is true, while the second is not. As presented in section 3.1, the shape of the free energy profile is not parabolic near the maximum (see Fig. 4). Therefore, we will use the more general formulation provided by eq. (15).

RESULTS

Free Energy

The free energy profiles for NP of various sizes were obtained using the methods described above. To gain a better conceptual understanding, the average bead density within the lipid bilayer was also calculated. Presented in Figure 4 is the bead density profile (a), accompanied by the free energy profiles for several of the studied NP (b). The free energy profile of the smallest NP, σ = 0.47, is qualitatively similar to the profiles obtained by fully atomistic simulations of small hydrophobic molecules, such as oxygen, preformed by Marrink.6

Figure 4 is the bead density profile (a), accompanied by the free energy profiles for several of the studied NP (b).

Figure 4. (a) Density profile of the DPPC bilayer. The solid line shows the average number of DPPC beads per cubic nanometer, while the dashed (dotted) line shows the average number of tail (head) group beads. (b) Free energy profile for three of the spherical NP studied. (c) Magnification of the circled region of (b). Observe the relatively low energy barrier for entering the bilayer.

As expected, all considered hydrophobic NP exhibit a preference for the bilayer interior over the bulk water phase. This preference is caused by the fact that hydrophobic NP have a greater attraction to the hydrophobic DPPC tails, as opposed to the hydrophilic head or water beads. As discussed in section 3.3, these stronger attractive interactions act as a driving force for the NP to remain in the bilayer interior. Also, we observe that larger NP exhibit a deeper free energy minimum, which is due to the greater opportunity for attractive interaction with lipid tails.

Figure 4c shows the details of the free energy profiles near the bilayer entrance. It is clear that relatively low energies, ~1 kJ/mol, are associated with entering the DPPC bilayer. Such a low energy barrier suggest that even larger particles, σ = 0.88 nm, have little resistance to entering the bilayer. Note that the complete shape of the free energy profile prevents the direct use of Kramers’s eq. (16) and requires the more general form, eq. (15), discussed in section 2.5.c.

As NP size increases, the development of local free energy maxima near the center of the bilayer is observed. The general region in which the local maxima reside is described by Marrink et al.6 as the soft polymer region. The region is characterized by the high density of ordered tail groups and was found to provide the most resistance for transport of small molecules. It is not surprising that larger particles encounter larger barriers in this region. More space for traveling is required for particles of larger size. In this region larger particles may find it difficult to push apart the ordered tail groups to provide sufficient space.

Diffusivity

The time-dependent friction coefficient at various positions in the bilayer was calculated using eq. (11) and (14). Figure 5 shows a few examples of the random force autocorrelation function for the smallest considered NP. From these results, it is clear that the correlation time is very small and, therefore, the Lagenvin equation with zero-correlated noise, eq. (8), is justified. Using eqs. (12) and (13), and the obtained γ, the nanoparticle diffusivity in the bilayer and is shown in Figure 6.

Figure 5. Examples of the autocorrelation function for the random force for the smallest considered NP

Figure 5. Examples of the autocorrelation function for the random force for the smallest considered NP (σ =0.47).

 

Figure 6. Diffusivity of the considered NP inside the lipid bilayer.

Figure 6. Diffusivity of the considered NP inside the lipid bilayer.

We note that the numerical integration of γ led to significant noise in the estimation of D. However, the average value obtained for the smallest NP, σ = 0.47 nm, is of the same order as the value found by fully atomistic simulations,6 for a particle of σ = 0.45 nm.

Figure 7 shows a decrease in diffusivity as NP size increases. This is to be expected, as larger particles commonly exhibit more friction than do smaller particles. All three NP exhibit a slight decrease in diffusion upon entering the bilayer. However, on the order of the noise the diffusivity inside the bilayer is relatively constant. This observation justifies the use of eq. (15) together with the assumption of constant diffusivity in our calculations of the NP transport rates.

Rate of Transport

We calculated the transport rate, <τtrn of NP through a bilayer using eq. (15). We took the point of entry, z1, to be ~ -3.5 nm and the point of exit, z2, to be ~ 3.5 nm (recall that z = 0 nm corresponds to the center of mass of the bilayer, see also Fig. 4). In this calculation we used the values of diffusivity averaged over all positions in the bilayer. The results are provided in Figure 7.

Figure 7. Semi-logarithmic plot of the mean transport time as a function of NP effective diameter.

Figure 7. Semi-logarithmic plot of the mean transport time as a function of NP effective diameter.

We observed that the average time to cross the bilayer increases significantly with NP size and that the increase is nearly exponential. This increase in transport time is caused by the increased friction, the deeper free energy minimum, and the added free energy barrier in the bilayer interior. Marrink et al.6 suggested that the bilayer interior acts as a trap for hydrophobic particles. Our study suggests further that the “trapping” effect is significantly increased with particle size.

Because the NP do not experience a significant energy barrier, we may assume they diffuse into the bilayer relatively quickly. This is confirmed by calculation of the mean escape time of a particle from the center of the bilayer to the bulk water phase. We computed the mean exit time using eq. (15) and found that it was approximately equal to the mean transport time. This shows that the NP enter the bilayer very quickly and spend the majority of the transport time inside the bilayer interior.

DISCUSSION

The results of this study yield a better understanding of the transport properties associated with carbon-based NP of various sizes through a DPPC bilayer. The free energy profile (Fig. 4) reveals that there is no significant energy barrier to enter the bilayer for any of the NP studied. In the absence of an energy barrier, the only resistance for our NP to enter a DPPC bilayer is the friction associated with diffusion. This would suggest, and is confirmed in Section 3.3, that these NP may enter the bilayer relatively quickly and that most of the transport time is spent inside the bilayer.

We observed that the NP are trapped by the hydrophobic interior of the bilayer. This effect is seen in the form of a free energy minimum (Fig. 4). Larger particles experience a lower free energy at this minimum, most likely due to an increase in the number of possible hydrophobic interactions. This suggests that size will play a crucial role in the time an NP will spend in the bilayer interior.

We also investigated even larger NP; however, the current model produced unreasonable results. As seen in Figure 8, the free energy profile of NP with σ = 1.5 nm exhibits a very large minimum, which corresponds to unrealistic transport rates. Such a result is most likely due to the limitations associated with modeling of such large NP by a single Lennard-Jones bead. Recall that we altered the value of the potential parameter σ to represent particles of various sizes. However, this change affects also the range of interactions. As σ increases, the range of attractive interactions may become unrealistic for van der Waals interactions, which would cause an unrealistically strong interaction of NP to lipid tails.

Figure 8. Free energy profile for a large NP, dotted line, compared to that of a smaller NP, solid line.

Figure 8. Free energy profile for a large NP, dotted line, compared to that of a smaller NP, solid line. The energy difference between bulk phase and bilayer interior is unreasonably large.

In the future, we plan to overcome this problem by modifying our model. Instead of using a single bead to represent large NP, it may be more realistic to use several smaller beads bonded together. This would allow us to increase the size of the NP without unreasonably increasing the range of attractive interactions. This alternative model would also allow us to explore the effect of molecular shape.

With the successfully modeled smaller NP, we were able to observe the size dependence on mean time transport time across the bilayer (Fig. 7). We observed that this average time increased nearly exponentially as a function of NP effective diameter. For the larger particle, σ = 0.88 nm, the mean transport time was relatively large, ~10-2 sec.

We might assume from these findings that larger NP (σ ~1.5 nm) will require even longer transport times. Because the typical diameters of fullerene molecules range from 1 to 1.5 nm, we might also infer that fullerenes will spend a significant time within the interior of nearby lipid bilayers, or cell membranes. This may allow more opportunity for such NP to impact the membrane interior, through physical interactions, or, if possible, chemical reactions. Therefore, it is possible that fullerene size plays a crucial role in their cytotoxic effect.

ACKNOWLEDGEMENTS

This work was supported by the Environmental Protection Agency, and by internal University of Florida Funds: SNRE and the University Scholars Program.

REFERENCES

  1. Venkatesan, N. ; Yoshimitsu, J.; Ito, Y.; Shibata, N.; Takada, K. Biomaterials 2005, 26, 7154.
  2. Oberdorster, E. Environmental Health Perspectives 2004, 4(10), 1881.
  3. Sayes, C. M.; Fortner, J. D.; Gue, W.; Lyon, D.; Boyd, A. M.; Ausman, K. D.; Tao, Y. J.; Sitharaman, B.; Wilson, L. J.; Hughes, J. B.; West J. L.; Colvin, V. L. Nano Letters 2004, 4(10), 1881.
  4. Finkelstein, A. J. Gen. Physiol. 1976, 68, 127
  5. Marrink, S. J.; Berendsen, H. J. C. J. Phys. Chem. 1994, 98, 4155.
  6. Marrink, S. J.; Berendsen, H. J. C. J. Phys. Chem. 1996, 100, 16729.
  7. Marrink, S. J.; Sok, R. M.; Berendsen, H. J. C. J. Chem. Phys. 1996, 104(22), 9090
  8. Smit, B. Phys. Rev. A 1988, 37, 3431
  9. Palmer, B. J.; Liu, J. Langmuir 1996, 12, 746
  10. Goetz, R.; Lipowsky, R. J. Chem. Phys. 1998, 108, 7397
  11. Shelley, J. C.; Shelley M. Y.; Reeder, R. C.; Bandyopadhyay, S.; Klein, M. L. J. Phys. Chem. B 2001, 105, 4464
  12. Marrink, S. J.; Mark, A. E. J Am. Chem. Soc. 2003, 125, 15233
  13. Marrink, S. J.; de Vries, A. H., Mark, A. E. J. Phys. Chem. 2004, 108, 750
  14. Goel, A.; Howard, J. B.; Vander Sande, J. B. Carbon 2004, 42, 1907
  15. van der Spoel, D.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Hess, B.; Feenstra, K. A.; Lindahl, E.; van Drunen, R.; Berendsen, H. J. C. Gromacs User Manual version 3.1.1, Nijenborgh 4, 9747 AG Groningen, The Netherlands. Internet: www.gromacs.org, 2002
  16. Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, NY, 1987
  17. Kopelevich, D. I.; Panagiotopoulos, A. Z.; Kevrekidis, I. G. J. Chem. Phys. 2005, 122, 044908
  18. Hummer, G.; Kevrekidis, I. G. J. Chem. Phys. 2003, 118, 10762
  19. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comp. Phys. 1977, 23, 327
  20. Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684
  21. Gardiner, C. W. Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Siences; Third Edition, Springer-Verlag: Berlin, 1983
  22. Zwanzig Nonequilibrium Statistical Mechanics; Oxford University Press: New York, NY, 2001
  23. Kramers, H. A. Physica 1940, 7 (4), 284

--top--

Back to the Journal of Undergraduate Research