Molecular dynamics simulations of carbon nanotube porins in lipid bilayers †

Artificial channels made of carbon nanotube (CNT) porins are promising candidates for applications in filtration and molecular delivery devices. Their symmetric shape and high mechanical, chemical, and thermal stability ensure well-defined transport properties, and at the same time make them ideal model systems for more complicated membrane protein pores. As the technology to produce and tune CNTs advances, simulations can aid in the design of customized membrane porins. Here we concentrate on CNTs embedded in lipid membranes. To derive design guidelines, we systematically studied the interaction of CNT porins with their surrounding lipids. For our simulations, we developed an AMBERand Lipid14-compatible parameterization scheme for CNTs with different chirality and with functional groups attached to their rim, and a flexible coarsegrained description for open-ended CNTs fitting to the MARTINI lipid model. We found that the interaction with lipid acyl chains is independent of the CNT chirality and the chemical details of functional groups at the CNT rims. The latter, however, are important for the interactions with lipid head groups, and for water permeability. The orientation and permeability of the pore are mainly determined by how well its hydrophobicity pattern matches the membrane. By identifying the factors that determine the structure both of isolated CNTs in lipid membranes and of CNT clusters, we set the foundation for a targeted design of CNT–membrane systems.


Introduction
The unique properties of carbon nanotubes (CNTs) have attracted a lot of interest, and a large number of possible CNT applications have been proposed. A particularly interesting property is their excellent water conductivity, as predicted from simulations 1 and conrmed in experiments. 2,3 Inserted as articial channels into lipid membranes [4][5][6][7] or assembled into densely packed hexagonal arrays, [8][9][10] CNT porins are candidates for drug delivery 11,12 and ltration. [13][14][15] Their transport 2 Computational methods

Atomistic simulations
We performed atomistic molecular dynamics simulations of single CNTs, pairs, and groups of four in lipid membranes. We created CNTs with four different functionalizations (COOH, COO À , OH and H-only) in zigzag and armchair congurations with three different lengths of the CNTs (4.0 nm, 4.5 nm 5.0 nm) and simulated them in a POPC membrane. To test for the effects of different lipids, we simulated a CNT of 4.5 nm with OH functionalization in membranes of DOPC, DPPC and POPE. In addition, we simulated pairs and quartets of CNTs in POPC using different starting conditions. We used the Lipid14 force eld 55 in Amber 14 (ref. 56) and derived GAFF-style 57-59 models for the differently functionalized CNTs.
2.1.1 Lipids and solvent. We modeled the lipids using the Lipid14 force eld. 55 Compared to previous attempts to model lipids with GAFF, 57 they were corrected such that simulations in a tensionless NPT ensemble give correct values for the area per lipid. 55,60 The starting congurations of the bilayers were generated using the soware Membrane Builder 61-63 via the graphical interface CHARMM-GUI. 64 With the help of the parameterization tools acpype 59 and antechamber, 58 generalized AMBER (GAFF) parameters can easily be created for other molecules, which are compatible with all AMBER force elds. The water molecules were modeled using the TIP3P model. About 150 mM NaCl was added. For all solvent atoms, the AMBER12SB force eld was applied.
2.1.2 Carbon nanotubes. The structures of the CNTs were generated using the python script BuildCstruct. 65,66 We extended this tool in order to include the possibility of adding the functional groups OH, COOH, or COO À at both ends of the nanotube. The modied code is available at https://github.com/bio-phys/cntgaff.
The GAFF parameters 57 were assigned to the nanotubes by running acpype 59 on these structures. Partial charges were dened by a standard scheme, which was parameterized as follows. For several nanotubes with different functionalizations and of different (small) sizes, the geometry was optimized and restraint electrostatic potential (RESP) charges 67 were calculated using R.E.D. (version III.52) with Gaussian 09 (RESP-A1: HF/6-31G* Connolly surface algo., 2 stage RESP t qwt ¼ 0.0005/0.001). The structures from BuildCstruct were either used directly as input, or (if the algorithm did not converge) a short atomistic energy minimization with non-optimized partial charges was performed.
We averaged and rounded the charges to two signicant digits (under the constraint to preserve the total charge of a functional group of 0 or À1, respectively) such that all identical functional groups have the same charge distribution and such that all charges beginning from the third carbon ring are assigned a value of zero. This is a good approximation according to our results as well as to results obtained by DFT calculations, 68 where partial charges of pristine single wall CNTs have been shown to concentrate on the edges and decrease rapidly towards the middle of the tube. Our results indicate that this is also true for functionalized CNTs.

Coarse-grained simulations
We simulated single CNT porins embedded in lipid membranes using Gromacs. 69 The lipid membrane was set up with the insane.py 70 tool. We varied the CNT in both length and diameter and repeated the study for three different functionalizations at the rim (SNda, SP3, and none), and simulated them all in a POPC membrane. In an additional study, we simulated one SNda-functionalized CNT in POPC/DOPC membranes of various size. Details on the procedure can be found in the ESI † and the force elds are described below.
2.2.1 The MARTINI lipid model. The coarse-grained (CG) simulations were performed using the MARTINI model. 71 In this model, in general, four heavy atoms are mapped to one CG interaction site (bead). For ring structures, there are special beads with a 2 : 1 mapping. Water is accordingly modeled with one bead for four water molecules.
Due to the attening of the energy landscape, the overall dynamics in coarsegrained simulations are accelerated. The developers of the MARTINI force eld recommend a scaling factor of four. 71 This means that the actual time is (approximately four times) longer than the simulated time.
2.2.2 Carbon nanotube model. We used a MARTINI-compatible CNT model generalized from the ones already used in previous work. 50,72 Our soware to generate the respective topology les and input structures for use with Gromacs is available online at https://github.com/bio-phys/cnt-martini.
For the carbon-only region of the nanotube, we chose the CNP interaction sites, which have already been used for fullerenes 73 and for capped nanotubes. 53 As in the latter model, end-group functionalization is modeled by replacing the CNP bead type with a more polar one. We used beads of the types SNda (weakly  polar) and SP3 (strongly polar). An example of this conguration is shown in Fig. 2. The bonded interactions between neighbor atoms were chosen as in the model for the capped nanotubes 53 with a bond length of 0.47 nm and a force constant of 5000 kJ mol À1 nm À2 .
For the angle potentials along the rings, we used a force constant of 350 kJ mol À1 rad À2 and an equilibrium angle of with N being the number of beads in one ring. We accounted for the stiffness along the tube by introducing harmonic improper dihedral angle potentials that maintained the angle b between two triangles with a force constant of 350 kJ mol À1 rad À2 , as shown in Fig. 2. As an equilibrium value, we get

Analysis
The python scripts used to analyze the simulations are publicly available. The scripts for the atomistic simulations can be found at https://github.com/bio-phys/ cnt-lipid14-analysis and those for the coarse-grained simulations at https:// github.com/bio-phys/cnt-martini-analysis. Both extensively use the MDAnalysis 74,75 python package. 2.3.1 Tilting. The tilting was measured using the main principle axis of the CNT. We calculated its angle with respect to the membrane normal, as given by the z axis of the system.
2.3.2 Two-dimensional radial distribution function. In analogy to the classical three-dimensional radial distribution function (RDF), we dene the twodimensional RDF g(r) using the number n(r, r + dr) of particles of a given type on average in a distance interval (r, r + dr) from an axis in space as gðr þ dr=2Þ ¼ nðr; r þ drÞ sA ring ðr; r þ drÞ with s being the area density of particles and A ring (r, r + dr) ¼ p(r + dr) 2 À pr 2 z 2prdr being the area of a ring around the center of the axis with radius r and thickness dr.

Order parameter.
The order of the lipids in a membrane is quantied using the deuterium order parameter dened as where q i is the angle between the bilayer normal and the vector joining the respective carbon atom C i to its H (resp. deuterium) atom, and h/i denotes an ensemble average. 55 For convenience, we report ÀS CD , averaged over all C-H bonds i in each acyl chain.
To estimate the inuence of the CNTs on the lipid ordering, we assigned the order parameter of each lipid acyl chain (averaged over all atoms) to the annular lipid shell to which the tail belongs. We used the minima of the RDF as cut-off distances to separate the shells. To get a bulk value for comparison, we simulated a patch of a pure lipid membrane of 288 POPC lipids. In coarse-grained simulations, a different denition of the order parameter has to be used. There are no hydrogen atoms and therefore no information on CH bonds. A rough estimate can be provided by the angle a j between the membrane normal and bond j connecting two MARTINI beads. This angle is supposed to be averaged over various C-C bonds and roughly perpendicular to the average orientation of a C-H bond. Due to the symmetry of cos 2 (x), the coarse-grained order parameter indicates lipid order in a similar way as ÀS CD . Again, we average over all bonds in each acyl chain. We note that this coarse-grained version and its atomistic counterpart cannot be compared quantitatively. However, we use them to compare general trends.

Water permeability.
The permeability of a channel can be estimated from uctuations in equilibrium simulations via a collective diffusion model. 76 This model was developed for CNTs of arbitrary diameter as a generalization of an earlier single-le model. 77 The model by Zhu et al. 76 monitors the water motion using the positions z i along the CNT axis of the S(t) water molecules in the CNT at time t. From the displacements dz i of each water molecule i˛S(t), the permeation change in an interval dt at time t is dened in a differential way as The net permeation n(t) is then obtained by integration with the initial condition n(0) ¼ 0. In equilibrium simulations, the trajectory of n(t) shows random-walk behavior. Its mean-squared displacement (MSD) is linear in time and denes a diffusion coefficient D n via An increase of n(t) by one means that one water molecule has been transported in the direction of positive z. Using the uctuation-dissipation theorem, 76 the pore permeability p f is obtained as 3 being the volume of a water molecule.

Interaction of single CNTs with lipids
The lipids form a strong annular shell around the CNT. We evaluated the radial distribution function (RDF) with respect to the principal axis of the CNT for the carbon atoms of the lipid acyl chains in all simulations containing only one CNT (Fig. 3). A pronounced ring-like multilayer structure is clearly visible, with a strong rst peak of the RDF followed by up to about ve further maxima. In accordance with recent simulations comparing different biological and biomimetic nanopores, 54 this structuring is a consequence, rst, of the tight lipid-CNT interactions and, second, of the relatively smooth surface of the CNT. Because we calculated the RDF with respect to the main CNT axis, to minimize the effect of CNT tilting, we could resolve the rst very pronounced peak in more detail. This peak shows a distinct shoulder on the far side. Inspections of the structures showed that this shoulder is caused by acyl chains adhering to the CNT surface only partly while the rest is kept close to the tube, but not in direct contact with it.
The amplitudes of all further maxima and minima decrease exponentially with the distance from the CNT (Fig. 3), and remain visible up to the h layer. In contrast to typical transmembrane proteins, the CNTs increase the lipid order in their proximity. For all CNTs, the order parameters in the rst four shells around them were signicantly higher than the respective bulk values. In the rst shell, they even became twice as high as in the bulk. Both effects, the pronounced lipid shells as well as the increase of the order parameter, were observed with minimal differences for the two CNT chiralities studied, and for every type of functional group at the rim.
In tests for possible artifacts of the nite box size (see ESI †), no dependence of the lipid layering on the box size was found. However, the increase in the order parameter near the CNT might be overestimated by z0.03-0.05 in the typical simulation boxes of this study (7-10 nm). In narrow boxes, lipid ordering around CNTs is more pronounced, likely enhanced by periodic boundary conditions. Therefore, caution is advised also in studies of lipid shells formed around membrane proteins in simulations using typical box sizes. Due to its regular hydrophobic surface and its mostly upright position, a CNT binds strongly to lipids and orders them in distinct rings. In the immediate vicinity of the CNTs, the lipids adopt structures not typical of a liquid-disordered phase. Instead, the lipids seem to form a strong coat around the nanotube, especially those in the very rst shell. Annular lipid shells are well-known for membrane proteins, 78,79 but the protein coating is not as strong and regular as seen here. 54 The effects of such lipid shells can be observed here in pure and unperturbed form. We can therefore consider CNTs as idealized models for the study of annular-lipid structures forming around membrane proteins (especially barrel-shaped proteins), and possibly as models for the induction of an "orderphilic" transition. 80 For different lipid types, the density of annular rings differs slightly but the main characteristics are preserved (Fig. 4). Those lipids that would be highly ordered even without the CNT show more pronounced maxima in the RDF (e.g., DPPC). The less-ordered lipids (e.g., DOPC) show a stronger shoulder in the rst peak. As this shoulder occurs for POPC but not for POPE, we conclude that it does not stem directly from an inherent property of the lipid acyl chains (which are the same in both lipids), but from the overall order. The most distinct behavior is seen for DPPC. Due to its two palmitoyl tails, it is substantially more ordered. It not only forms a dense layer around the CNT but also all further density peaks are sharper and their distance is a bit shorter than for the other lipids. By and large, however, the behavior is qualitatively the same for all lipid types.

CNT functional end groups
In contrast to the lipid acyl chains, the RDFs of the lipid head groups show a strong dependence on the rim functionalization of the CNT. For armchair and zigzag CNTs, the RDFs show a common trend (Fig. 5). The head-group density at the edge of the CNT increases with the length of the CNT. The longer tubes order the polar head groups around them more strongly, and keep them away from the CNT openings. For short tubes without or with only weakly polar end-  functionalization, the head groups can bend over the edge of the CNT (Fig. 6), especially when the CNT is tilted. This effect is represented in the RDF as a peak in the region where r is smaller than the radius of the CNT. CNTs with  carboxylate (COO À ) functionalization, carrying a negative charge, manage best to prevent the head groups from bending over the rim and into the opening. Depending on the purpose of the membrane, charged rims might be more or less desirable. In any case, polar functional groups can prevent the pore from closing by lipid blockage.

Comparison to coarse-grained simulations
The structure of the lipid shells seen in the atomistic simulations is well reproduced in the coarse-grained MARTINI model simulations (Fig. 7). The distances and widths of the distinct shells are in excellent agreement. Lipid structuring extends out to ve layers in both models. Even though the rst peak is less sharp in MARTINI (which is to be expected in a coarse-grained model), it is still higher than that extrapolated from the exponential behavior of the outer shells. The MARTINI model is therefore suitable to model the effects of CNT-lipid-tail interactions.
The order parameters in the MARTINI model show the same trend as in the atomistic model, but it is less pronounced. The order parameter is strongly increased in the rst two shells around the CNT and slightly increased up to the h shell. The difference between the models is not dramatic as, by denition, the MARTINI order parameter is not able to capture the ordering effect on single CH bonds, but only the much coarser general orientation of the tail.
As in the atomistic simulations, no dependence of lipid layering on the box size was found for the MARTINI model (see ESI †). In contrast to the atomistic model, no trend was found for the order parameter. This is not surprising as the overall ordering effect is less pronounced in the MARTINI model. The box sizes in the MARTINI simulations (z15 nm) are therefore sufficient.
We note that the CNT chirality and the specic details of the functional groups are not resolved within our coarse-grained MARTINI model. However, the  MARTINI CNTs reproduce well the structural properties found in the atomistic simulations. These are dominated by geometry and the hydrophobic effect.

Tilting of CNTs
As a general trend in the atomistic simulations, the longer tubes tend to tilt more with respect to the membrane normal than the shorter ones and polar functional groups prevent tilting (Fig. 8). These two observations can be explained by the CNTs striving to minimize hydrophobic mismatch, i.e., the contact of the hydrophobic CNT surface with water or hydrophilic lipid head groups. Even though we can identify trends within a large number of simulations, the individual simulations are too short for tilting to be fully equilibrated. A study on one CNT (armchair, 5 nm, COOH) in differently sized boxes (see ESI †) indicates that the tilting is only marginally affected by the nite size and rather under-than overestimated in small boxes. To achieve much better sampling in fully equilibrated simulations, we turned to coarse-grained methods. To study the geometry effects on tilting (Fig. 9), we compared CNTs of three different diameters and ve different lengths, each for one type of functional end group: strongly polar (SP3), weakly polar (SNda), and no functionalization at all. We conclude from the Fig. 8 Tilting angle of the CNTs with respect to the z axis in atomistic simulations of single CNT porins in a POPC lipid membrane. The chirality (armchair, zigzag) was varied as well as the length of the CNT (4.0 nm, 4.5 nm, and 5.0 nm). Rims with four different functionalizations (COO À , COOH, OH, none) were studied. Vertical bars denote the standard deviation. MARTINI box-size study (see ESI †) that for the box sizes studied we obtained an unperturbed estimate of tilting. The coarse-grained tubes conrm the trend from the atomistic simulations: long tubes tilt, especially when they are thin and have no polar groups at their rim. Below 4.5 nm length (the thickness of the membrane), the CNT length has almost no inuence. These results are in full accordance with the atomistic ones and conrm that hydrophobic mismatch is the main determinant of the CNT porin tilt.

Water structure and permeation
Water in nanopores is strongly conned. In our simulations, the water inside the CNTs adopted cylindrical shapes (Fig. 10), which are more pronounced in the wider armchair CNTs. The radius of the volume accessible for the water (carbon   coordinates minus van der Waals radius of a carbon atom) is 3.8Å in our armchair CNTs and 3.0Å in our zigzag CNTs. The functional groups attached to the rim can signicantly change the permeability of a CNT porin (Fig. 11). Whereas for the wider armchair CNTs in this work the permeability varies considerably, a clear trend is visible for the narrower zigzag CNTs. The polar hydroxyl (OH) and carboxyl (COOH) groups clearly enhance permeability compared to an unfunctionalized CNT. They prevent tilting and blockage of the pore by lipid headgroups, as discussed in the previous sections. The negatively charged COO À groups have the same effects structurally, but still reduce water permeability. They attract cations (here Na + ) that reside at the pore entrance and partially block it. In the smaller CNTs (here those with a zigzag conguration), this leads to an almost complete blockage of the pore. These results show that, especially for small radii, rim functionalization is an important factor determining the permeability.

Carbon nanotube pairs and quartets
In atomistic simulations of more than one nanotube, lipids formed strong layers around the tubes for CNT pairs as well as quartets. Order is especially strong between two CNTs (Fig. 12). The lipid acyl chains stretch out and adapt to their squeezed-in position. This results in structures that are quite stable at least on the time scales of our simulations. The resulting conguration of CNTs in each simulation strongly depends on the starting conditions (see the ESI † for more details). CNT pairs that were prepared in direct carbon-carbon contact stayed in contact and pairs prepared with a distance of 2 nm between the central axes of the CNTs remained separated by one layer of lipid. Quartets of CNTs that were initially prepared in a square conguration with 2.5 nm distance remained separate. Quartet bundles initially prepared in contact quickly turned towards a stable rhombic shape (Fig. 13,  right). This represents the hexagonal lattice structure expected on large scales due to the cylindrical shape. 8 The lipid order parameter in all structures with more than one CNT (Fig. 14) shows the same behavior as for single CNTs. In our analysis, we distinguished shared (S) lipids that belong to the rst shell of two CNTs and are therefore in a particularly conned conguration. For them, the order parameter varied more strongly due to the lower number of samples and depending on the exact position in which the lipid got trapped. However, the trend over all simulations shows a slight increase in order. Overall, we expect a lipid membrane densely packed with CNTs to be more ordered and maybe to undergo a phase transition to a gel phase.

Conclusions
We performed atomistic simulations of single CNT porins with different lengths, functionalizations and chiralities, and of CNT pairs and quartets. We also performed simulations of membrane embedded CNTs using a coarse-grained MARTINI-compatible model. In both atomistic and coarse-grained models, we observed pronounced annular lipid shells around the CNTs. The inuences of  functionalization and CNT chirality on the shell are negligible. The strong structuring of the lipids is therefore a universal feature for CNT porins. We observed in our simulations that the position and stability of the porin within the lipid bilayer depended on its length and the hydrophobicity of the functional groups at the rim. We found in atomistic and in coarse-grained simulations that CNTs minimize water exposure of their hydrophobic surface. In recent experiments, 5 CNTs were observed without signicant tilt. The mechanism behind this stability is still unclear. The tendency to tilt might be counterbalanced by external factors that minimize the exposure of the hydrophobic CNT surface, such as coating lipids or detergents, or by mechanical effects not captured in our simulation setup.
To control the openings of a CNT porin, different kinds of polar rim functionalizations can be used. Polar groups help prevent the blocking of short porins by lipid head groups, keep the pore open, and thus increase the efficiency of CNTmediated transport across the membrane. Charged groups, however, attract ions that can block narrow pores and lower the water permeability. As a consequence, polar end functionalization affects the embedding structure and permeability of lipid-membrane embedded nanopores and in turn modulates the transport properties, e.g., by making the pore selective. 21 The formation of CNT bundles within the membrane or the pairing of CNT porins result in stable membrane-spanning congurations with pores open to the solvent. Packing of CNTs into lattice structures stabilizes their vertical orientation. The lipid layering effect around these axially aligned CNTs is at least as strong as around isolated porins, and it might even have a stabilizing and structuring effect on the assemblies. We thus expect a large concentration of CNT porins to increase the permeability per CNT. However, such CNT bundles or more extended arrays will alter the mechanical properties of the membrane.
In conclusion, we demonstrated that CNT porins exert a strong inuence on the surrounding membrane lipids. General principles emerge that articial water channels made from CNT porins should fulll, providing guidance for further experiments.

Conflicts of interest
There are no conicts to declare.