Elucidating collective translocation of nanoparticles across the skin lipid matrix: a molecular dynamics study

Abstract

The top layer of skin, the stratum corneum, provides a formidable barrier to the skin. Nanoparticles are utilized and further explored for personal and health care applications related to the skin. In the past few years, several researchers have studied the translocation and permeation of nanoparticles of various shapes, sizes, and surface chemistry through cell membranes. Most of these studies focused on a single nanoparticle and a simple bilayer system, whereas skin has a highly complex lipid membrane architecture. Moreover, it is highly unlikely that a nanoparticle formulation applied on the skin will not have multiple nanoparticle–nanoparticle and skin-nanoparticle interactions. In this study, we have utilized coarse-grained MARTINI molecular dynamics simulations to assess the interactions of two types (bare and dodecane-thiol coated) of nanoparticles with two models (single bilayer and double bilayer) of skin lipid membranes. The nanoparticles were found to be partitioned from the water layer to the lipid membrane as an individual entity as well as in the cluster form. It was discovered that each nanoparticle reached the interior of both single bilayer and double bilayer membranes irrespective of the nanoparticle type and concentration, though coated particles were observed to efficiently traverse across the bilayer when compared with bare particles. The coated nanoparticles also created a single large cluster inside the membrane, whereas the bare nanoparticles were found in small clusters. Both the nanoparticles exhibited preferential interactions with cholesterol molecules present in the lipid membrane as compared to other lipid components of the membrane. We have also observed that the single membrane model exhibited unrealistic instability at moderate to higher concentrations of nanoparticles, and hence for translocation study, a minimum double bilayer model should be employed.

---

1 Introduction

The outermost layer of the skin, stratum corneum (SC), provides a barrier against chemical, biological and mechanical stimuli.^1,2 This layer allows only small and relatively lipophilic molecules to diffuse into the deeper layer of the skin.³ The thin layer (∼20 μm) is composed of corneocytes and a lipid matrix, where the latter one fills the extracellular space, and the whole assembly resembles a brick-and-mortar structure.⁴ The passive permeation of an external entity (molecules, particles, and formulations) is mainly controlled by the lipid matrix, and it is indeed responsible for the barrier function of the skin.^3,5 A disturbance in the lipid composition and its lateral packing can lead to disease conditions such as atopic dermatitis and psoriasis.⁶

A healthy skin barrier is necessary for defending the skin from the attack of foreign entities (pathogens, bacteria, and harmful substances), but at the same time it creates a hurdle in delivering therapeutics on/through the skin.^7,8 Hence, several strategies are being developed and used to breach the skin barrier selectively and safely for numerous healthcare and personal care applications.^9–11 One of the strategies is to apply an external energy source (electric field,¹² ion flux,¹³ and ultrasound¹⁴) or mechanical source such as microneedles¹⁵ to name a few. These methods require not only dedicated devices but also sophisticated control of operating parameters (electric field, voltages, pulse rate, microneedle size, etc.). The other strategies focus on modifying the therapeutic formulation itself by the incorporation of permeation enhancers,¹⁶ skin penetrating peptides,¹⁷ and nanomaterials.¹⁸

The permeation enhancers are a class of molecules known to enhance the skin permeation by either extraction of lipids or fluidization of the skin matrix.^19–21 The penetration enhancers not only affect the lipid matrix but also the skin proteins, which could lead to irritation and may have toxic effects.^19,20 Ethanol and dimethyl sulfide oxide are two common permeation enhancers, and it has been shown both by experiments^22,23 and simulations^24,25 that they can disrupt the lipid packing to undesirable limits at higher concentrations.

Researchers have also studied nanoparticles (NPs)/nanomaterials as they have unique chemical and physical size dependent properties and can be tailor made from organic (lipids and biopolymers) and inorganic materials (metals).²⁶ Based on the application one can choose either deformable nanomaterials like liposomes or rigid metallic nanoparticles. Liposomes are believed to cross the skin barrier by its deformation and merging with the skin lipid membrane; however exact mechanisms are still unknown.²⁷ Gold NPs have been around for the past few decades and been explored for various biological applications due to ease of modification of their surface with chemical or bioactive molecules to ensure the binding of different types of substrates over their surface.²⁸ Gold NPs of various shapes, sizes and surface morphologies have been synthesized and shown to penetrate the deeper layer of the skin.^29–33 Not only experiments but also molecular simulation studies have been carried out to understand the permeation mechanism of NPs across the cell membrane^34,35 and skin lipid layer.^36,37 The in silico models and methodology have been utilized for the design of NPs of various morphologies to overcome the barrier function of skin and the intestinal membrane^38,39 as well. The potential of mean force or permeation free energy of single nanoparticles across the skin lipid bilayer model was calculated using umbrella sampling simulations. The free energy of permeation of bare and coated nanoparticles for a single bilayer system exhibited a minimum in the middle of the bilayer system. The coated and bare nanoparticles were observed in the middle of the bilayer in a long unrestrained MD simulation run as well. In the skin SC, multiple lipid layers are present and in experiments it has been shown that nanoparticles do cross these layers. Based on single nanoparticle permeation free energy and unrestrained simulation these phenomena are difficult to explain.

The simulation studies, so far presented in the literature, focus on the permeation of single nanoparticles across either the skin lipid or the cell plasma bilayer.^34–37 However, several NPs are present in a formulation and when applied on skin, it is highly unlikely that nanoparticle-skin lipid interactions are not being affected by the other nanoparticle–nanoparticle interactions. There are some studies which have reported the co-operative transport of NPs across the cell membrane.^40,41 Yue and Xianren Zhang reported the cooperative effect of nanoparticle permeation across the model lipid bilayer using dissipative particle dynamics (DPD) simulations.⁴⁰ It was shown that small NPs formed small clusters on the membrane and internalized as a whole entity, whereas the large nanoparticles remained dispersed and internalize independently. Four types of internalization pathways, namely, synchronous internalization, asynchronous internalization, pinocytosis-like internalization, and independent internalization, were observed in simulations and found to be dependent upon the NP size and their inter-particle distances. Zhang et al. performed DPD simulations of a model membrane with NPs of various surface chemistry/patterns.⁴¹ The translocation of NPs across the bilayer was observed in a cooperative manner and highly influenced by the surface chemistry of the beads of NPs. The hydrophilic NP permeation increased with their concentration, but for partial and totally hydrophobic NPs, the opposite trend was reported.

Unlike the cell membrane, the skin lipid matrix is multilayer in nature and the permeation mechanism could be very different from that of a single cell membrane. Moreover, in the above studies,^40,41 the receptor mediated translocation process was modelled by modifying the interaction parameters, whereas the transport across the skin SC layer is not receptor mediated. Hence, passive permeation of NPs of various sizes and surface chemistry across a multilayer model membrane would be entirely different. In this study, we report the translocation of NPs of various sizes and surface chemistry across the multilayer skin lipid matrix model. Unrestrained coarse grained (CG) molecular dynamics (MD) simulations were performed for bare and thiol coated NPs. Two kinds of bilayer systems, namely single lipid bilayer (SLB) membrane and double lipid bilayer (DLB) membrane, were simulated, and it is shown that the DLB membrane model should be used to study the translocation of nanoparticles across the skin lipid membrane.

2 Methods, systems and models

The permeation of molecules through skin lipid layers generally occurs on millisecond (ms) to micro-second (μs) time scales depending upon the size, shape, and surface chemistry of molecules. To capture the permeation events on realistic time and length scales, CG models derived from the MARTINI force field were used.^42–44 The NP parameters were taken from earlier studies.^45,46 The gold nanoparticle was taken out from an FCC crystal of a nanoparticle. A one-to-one mapping from an atomistic atom to a CG bead was used. The NP beads were modelled as “C5” MARTINI bead type. The skin SC lipid lamella is made up of three classes of lipids namely ceramides, cholesterol, and free fatty acids. Ceramides can be categorized into 18 subclasses based on the number and positions of hydroxyl groups at the lipid head groups.^2,6 The human skin SC typically has an equimolar ratio of ceramides (CER), free fatty acids (FFA), and cholesterol (CHO).¹ The CG structure of the skin lipid bilayer (equilibrated for 3 μs) was taken from earlier work.³⁹ The bilayer system had 936 each CER, CHO and FFA and 11 520 CG water beads. The bilayer size was 19.88 nm × 19.88 nm × 11.87 nm. Here onwards this bilayer system is termed the single lipid bilayer (SLB) membrane model. Another model termed the double lipid bilayer (DLB) membrane model was prepared using the SLB. The water layer from the equilibrated SLB was removed, and the remaining lipids were replicated in the z direction and solvated again with water at both ends. The system was subjected to energy minimization and subsequently NVT and NPT runs.³⁹

To allow complete inclusion of nanoparticles in the upper water layer and to avoid artifacts associated with the system size and nanoparticle–nanoparticle interactions over periodic boundaries, the water content was increased on each side of the bilayer. The NPs were inserted in the upper part of the water layer at a distance of ∼6–7 nm from the center of the lipid bilayer. The overlapped water molecules were removed, and systems were energy minimized. The systems were further equilibrated for 200 ns each under NVT and NPT conditions by keeping lipids and NPs fixed to ensure proper solvation. Properly solvated and equilibrated systems were then run for 1 μs in an NPT ensemble.

All simulations were carried out in NVT and NPT ensembles using the GROMACS MD package.^47–50 The temperature was controlled at a skin temperature of ∼305 K, using the Berendsen (equilibration run) and Nose–Hover (production run) thermostats with a time constant of 2 ps. Pressure was controlled by using Berendsen (equilibration run) and Parrinello–Rahman (production run) barostats with a time constant of 6 and 12 ps respectively and compressibility of 4.5 × 10⁻⁵ bar⁻¹ with semi-isotropic coupling. The pressure was independently controlled in XY and Z directions to obtain a tensionless lipid bilayer. The LJ potentials were smoothly shifted to zero between a distance rshift = 0.9 nm and a cutoff distance of 1.2 nm. The pair list was updated at every 20 steps. The configuration was sampled every 100 ps in the production run.

3 Results and discussion

3.1 Single lipid bilayer membrane model

The SLB models have been used to study the translocation of the nanoparticle across the cell and skin membranes.^34–37 To check whether the same model can provide insights into the permeation or translocation of multiple NPs, simulations with two types of NPs were performed. The bare nanoparticle had a surface bead type of “C5” of MARTINI force field. The bare NPs coated with dodecane thiol were termed coated NPs. The details of the systems simulated are given in Table 1. The systems were labeled ‘B’ and ‘C’ for bare and coated NPs respectively.

Table 1 Details of the simulated systems

System	Membrane	Number of NPs	Size (nm)	Type	Simulation time(μs)
a Extended runs.
B2	SLM	2	3	Bare	1
B3	SLM	3	3	Bare	1
B4	SLM	4	3	Bare	1
B5	SLM	5	3	Bare	1
B6	SLM	6	3	Bare	1
B7	SLM	7	3	Bare	1
B8	SLM	8	3	Bare	1
B9	SLM	9	3	Bare	1
C2	SLM	2	3	Coated	1
C3	SLM	3	3	Coated	1
C4	SLM	4	3	Coated	1
C5	SLM	5	3	Coated	1
C6	SLM	6	3	Coated	1
C7	SLM	7	3	Coated	1
C8	SLM	8	3	Coated	1
C9	SLM	9	3	Coated	1
B13	DLM	3	1	Bare	1 + 3^a
B16	DLM	6	1	Bare	1 + 3^a
B19	DLM	9	1	Bare	1 + 3^a
B33	DLM	3	3	Bare	1 + 3^a
B36	DLM	6	3	Bare	1 + 3^a
C13	DLM	3	1	Coated	1 + 3^a
C16	DLM	6	1	Coated	1 + 3^a
C19	DLM	9	1	Coated	1 + 3^a
C33	DLM	3	3	Coated	1 + 3^a
C36	DLM	6	3	Coated	1 + 3^a

---

NPs of diameter 3 nm were utilized for studying the interaction with the SLB model. MD simulations were performed by initially placing the NPs at the top of the SC membrane (Fig. 1A and 1B) as described in Section 2. As summarized in Table 1, simulations were performed for 1 μs with 2, 3, 4, 5, 6, 7, 8, and 9 NPs. The snapshots of the systems at the end of the simulation run are shown in Fig. S1 (see the ESI†). Spontaneous partitioning of NPs was observed in all simulations (Fig. 1, S1 and S2†), irrespective of the surface coating and concentrations. NPs have exhibited events of partitioning as individual monomers (Fig. 1A and S1†) as well as for the aggregates (Fig. 1B and S1†). These findings highlight a highly favorable environment offered by the SC membrane for the NPs. Membrane configurations indicate alteration in the membrane packing after the partitioning of the NPs, and the distortion appears to be more dominating for the coated NPs (Fig. S1†).

Fig. 1 Bare and coated nanoparticles partitioning into the SLB model. Side-view MD snapshots of the (A) B2 and (B) B9 systems having initial (top) and final configurations with and without lipid molecules. Time trace for the z-coordinates of the center of mass of NPs, and CER headgroups of upper and lower leaflets for (C) B2 and (D) B9. The CER, FFA, CHO and NPs are represented in red, blue, green, and yellow colors respectively. The images were rendered using VMD software.⁵¹

The above results are in line with some of the experimental evidence of finding the gold NPs in the deeper layer of the skin.^29–33 Huang et al. showed that gold NPs could deliver the model protein drugs in the skin dermis.²⁹ The findings of hydrophobic NPs in the interior of the lipid bilayer have been reported in some of the earlier simulation studies of plasma membranes as well.^35,37 Fullerene like NPs (size ∼1.2 nm) also have been found in the interior of the phospholipid bilayer and it's been shown that the NPs remained dispersed in the interior of the bilayer even at moderate concentrations of NPs.^34,52–55 The bare (3 nm) and coated (∼5 nm) nanoparticles are bigger than the bilayer width and found to be agglomerated even at moderate concentrations (Fig. S1†). This implies that not only the concentration of NPs but also the NP size determines their organization in the membrane environment.

We have analyzed the trajectories of the center of mass of ceramide headgroups (upper and lower) and the center of mass of each NP along the z-direction to investigate the partitioning behavior of the NPs (Fig. 1C, D and S2†). For the B2 system with two bare NPs (Fig. 1A), one NP was initially placed at the membrane–water interface while the other was placed above the SC membrane and in bulk water (Fig. 1A and C). Spontaneous insertion was observed for both the particles (Fig. 1A and C), while the first particle got partitioned within the initial 100 ns, the other one in bulk water entered the membrane after 300 ns. Once inserted, the NPs were present in the membrane core for the rest of the trajectory. The swelling of the membrane in the vicinity of NPs is quite apparent from the configurations (Fig. 1A) due to the hydrophobic mismatch between the membrane thickness and the size of the particles.

For the B9 system, the higher concentration of the NPs resulted in the aggregation in bulk water (Fig. 1B). Interestingly, the membrane allows the partitioning of this aggregate and accommodates a high concentration of NPs (Fig. 1B). Although membrane packing appears to be significantly distorted due to the lipid reorganization in the vicinity of the NPs, the overall bilayer structure is sustained. From the trajectory time trace data (Fig. 1D and S2†), it can be observed that five NPs were part of the aggregate, which is evident from the strongly correlated z-coordinates of those particles and insertion from the upper leaflet. However, one NP got partitioned into the membrane within the initial 100 ns from the downside. Two NPs formed another aggregate that entered the membrane environment from the lower leaflet. Although both aggregates entered the membrane at quite close time points, the membrane accommodated all particles while maintaining the bilayer structure (Fig. 1B).

All other systems (except B3 and B8) with bare and coated NPs exhibited similar behavior (Fig. S1†) where the NPs partitioned into the SC membrane within 500 ns (Fig. S2†). The B3 system presented an exception where the particles formed an aggregate in the initial stage and got partitioned into the membrane towards the end of the 1 μs simulation (Fig. S2†). Nevertheless, the partition was witnessed in all the simulations. In B8, the distortion is quite significant, as apparent from the membrane configuration (Fig. S1†). Higher distortion in the membrane topology in the presence of coated particles compared to the bare particles (Fig. S1†) is evident from the time trace data for headgroups, which show significant displacement in the membrane trajectories for coated particles.

The bilayer structure of the SC membrane has been altered due to the interactions with the NPs. The reassembly of lipids and distortion in the membrane is apparent from the membrane configuration corresponding to the last frame of 1 μs simulations (Fig. 1B and S1†). To quantify the perturbation in the membrane system due to the presence of a larger number of NPs, the head group fluctuations were evaluated. We have calculated the distributions of the CER headgroup fluctuations for the entire 1 μs trajectory and the last 100 ns of the trajectory to comprehend the alteration in the membrane topology (Fig. 2 and S3†). It is evident from the headgroup distribution for 1 μs data that the coated NPs have significantly perturbed membrane headgroups (Fig. S1†). The membrane has shown significant displacement for C5 as compared to the bare counterpart B5 (Fig. 2A and S2†). Also, the spread for the coated particles (∼5 nm) is significantly larger than the bare particles (∼3 nm). Higher fluctuations in the presence of coated particles compared to the bare particles are intuitive due to the presence of protruding thiol groups in the former.

Fig. 2 SLB model exhibits unrealistic headgroup fluctuations. Distribution for CER headgroup fluctuations for B5 (gold) and C5 (black) for the (A) entire trajectory and (B) last 100 ns. (C) Distribution for CER headgroup fluctuations for B13 (gold) and C13 (black) for the entire trajectory (top) and last 100 ns (bottom).

Final configurations suggest a stable membrane for both bare and coated NPs (Fig. S1†). Consequently, we have also determined fluctuations for the last 100 ns trajectory (Fig. 2B and S3†). Interestingly, it can be perceived that the fluctuations are settled towards the end of simulations for both cases, highlighting a tendency of the membrane to accommodate both types of NPs. The fluctuations are of the same order for both bare and coated NPs. During insertion, it can be concluded that perturbation in membrane headgroups for coated NPs is significant compared with the bare NPs (Fig. 2A and S3†). Although once the membrane gets equilibrated after accommodating the particles, the fluctuations become similar for both particles (Fig. 2B and S3†). These trends of higher distortion in the presence of coated NPs are conserved for all the cases (Fig. S1 and S3†), except for B8, where the membrane is entirely distorted (Fig. S1†), and hence a head-on comparison is not feasible. Also, all the simulations manifest undulated membranes towards the end of trajectories (last 100 ns).

The presence of NPs cause disturbances in the lipid packing, as can be seen from the projected area of the bilayer on the XY plane (Fig. S4†). Irrespective of the nanoparticle type (bare or coated) the area decreased with increase in the number of NPs; however for any given number of NPs, the projected area was smaller for the coated nanoparticle as compared to the bare one. This implies that coated NPs induced more undulation in the bilayer as compared to the bare NPs. The stability of the bilayer is also evident from the density profiles shown in Fig. S5.† The presence of dual peaks in the lipid's density corresponds to the existence of the bilayer structure. The dual peaks disappeared with an increased number of NPs for both bare and coated cases. The swelling of the bilayer is also apparent from the increased spread of density in the z direction. The tail order parameter (Fig. S6†) again confirms that the presence of a larger number of NPs significantly affected the lipid packing. The order parameter decreased with the increased nanoparticle number and coated nanoparticles had a much pronounced effect on the order parameter as compared to the bare nanoparticles for a given concentration of NPs. The bare NPs created small clusters in the interior of the membrane, whereas coated NPs created a single large cluster (Fig. S7†-Fig. S10†). Hence, membranes with coated NPs were found to be in a highly undulated form.

The SC membrane comprises three distinct types of lipids: CER, FFA, and CHO. To understand if the translocation of the NPs with the SC membrane is being stimulated by any selective lipid or all of them, we have examined their interactions independently. We have perceived all three lipids in the proximity of the NPs (Fig. 3A), which highlights strong interactions of the NPs with all three lipids. To understand the correlation within all components, we have evaluated the lateral density of all four components: NPs, CER, CHO, and FFA. The high lateral density of CHO molecules in the vicinity of the NPs exhibits a strong tendency of CHO molecules to bind with the NPs (Fig. 3B and S7–S10†). Although this correlation contrasts in the case of the CHO molecule it also persists for CER and FFA (Fig. S7–S10†). Interestingly these trends appear to be consistent for both types of particles (Fig. S7–S10†).

Fig. 3 Lipid molecules show selective binding with NPs. (A) MD snapshot of the lipid molecule within 1 nm of a gold NP. Lateral density map of (B) NPs and (C) CHO for B9. (D) Normalized contacts of lipid species with the NPs. Contacts are calculated for the last 100 ns of trajectories and normalized with the number of MARTINI beads corresponding to each lipid species and the number of NPs in the system. The CER, FFA and CHO are represented with red, blue, and green colors, respectively.

The aggregation of CHO around the nanomaterial has been reported before.^56–58 Guo et al.⁵⁶ reported the destabilizing effect of CHO on the lipid membrane due to the strong CHO–CNT interaction, and cholesterol triggers the translocation of an ultrashort CNT through the cell membrane. Raczyński et al.⁵⁷ showed strong interaction between CHO and silicon-carbide nanotubes, due to which the nanotube could be used for effective extraction of CHO localized on a surface of the endothelial protein 1LQV. In another study by Raczyński et al.,⁵⁸ they showed differences in CHO dynamics near nanotubes based on whether they are made of carbon or silicon carbide in the presence/absence of water.

To quantify these interactions and have a more robust gauge for comparison between lipid species, we have performed contact analysis over the last 100 ns of trajectory. We have calculated total contacts made from NPs with the lipid species, where contact is defined when a lipid CG bead is present within the radius of 1 nm of an NP. The normalized contacts (Fig. 3D) are calculated by normalizing the total contacts over the trajectory by the number of NPs in each system and the number of beads present in individual lipid species. A strong interaction of NPs with all three components can be recognized, and hence the NPs partition into the membrane. We observe a preferential binding of CHO molecules compared to CER and FFA in all the cases (Fig. 3D). The strong interaction of NPs with CHO highlights the capability of CHO as an effective receptor for NPs. For systems corresponding to bare NPs, the differences between CHO and the other two species are remarkably distinct. However B3 shows a considerably lower contact due to the slow permeation process than the other cases (Fig. S2†). It can also be observed that the bare NPs have a stronger tendency to interact with CER when compared with the FFA. The trends for normalized contacts are quite reproducible across all the cases showing a strong preference of CHO binding with the NPs irrespective of concentration and aggregation.

For coated NPs, there is a notable difference in the trends of contact analysis. The total contacts are on the higher end for coated particles than the bare particles due to the protruding thiol groups, which provide additional interaction sites. Particularly relative intensification in FFA contacts can be recognized. The tendency of preferential interactions of bare NPs with the CER compared to the FFA gets interchanged in the case of coated NPs, where the FFA shows an enhanced binding in the vicinity of the NPs. Indeed, CER and FFA contacts are comparable to the CHO contacts for coated NPs, while there was a significant difference in the case of bare NPs. The finding of the NPs in the vicinity of the CHO can also be seen from the radial distribution function (Fig. S11†) plots.

In Section 3.1, the results of interaction of NPs of various size and surface chemistry with the SLB model at various concentrations of NPs are discussed. A high concentration of NPs induced several structural changes in the bilayer and one reason could be the presence of a water layer at both ends of the bilayer. In the skin SC, several lipid layers are arranged in the lamellar form.

A SLB model for the SC membrane is deficient in representing the multilayer architecture of this membrane. To ascertain whether the trends addressed hitherto are realistic and not an artifact due to the simplistic SLB model, we have employed a multilayer (two layers) model to represent the skin barrier, which is closer to real systems. Also, to understand the effects of NP size on the partitioning properties, we have studied two sizes of NPs (1 nm and 3 nm) using this DLB (Fig. 4A and B). We have studied three systems for 1 nm particles with three, six, and nine NPs (Table 1 and Fig. S12†). For 3 nm particles, two simulations were performed using three and six NPs (Table 1 and Fig. S12†). Two bilayer simulations were carried out in a stepwise manner, where in the first step NP partitioning was observed for a duration of 1 μs. The water molecules were reduced in the second step to reduce the system size, and simulations were extended for another 3 μs. The assertive partitioning behavior observed in the SLB model was effectively reproduced in the multilayer model (Fig. 4 and S12†). Like what was discerned in the SLB model, the NPs were partitioned into the DLB SC membrane within the initial 600 ns of simulations, except for C33 (Fig. S13 and S14†).

Fig. 4 Bare and coated nanoparticles partition into the DLB model. Side-view MD snapshots of the (A) B13 and (B) C13 systems having NPs and a double-bilayer SC model corresponding to the initial (top) and final (bottom) configurations. Lipids are hidden in some cases for clarity. Time trace for the z-coordinates of the centre of mass of NPs, and CER headgroups of the upper, middle, and lower leaflets for (C) B13 and (D) C13. The same colour scheme is used for both plots. The CER, FFA, CHO and NPs are represented in red, blue, green, and yellow colours respectively. The NPs are highlighted in circles. The images were rendered using VMD software.⁵¹

In B13, two NPs were initially placed in the headgroup region of the membrane, and the third one was placed in bulk water above the DLB model (Fig. 4A). All three particles were restricted to the headgroup region of the membranes for the entire trajectory (Fig. 4A and C). It can be concluded from the trajectories of the 1 nm bare particles in all simulations that a favorable site is present for such molecules in the headgroup region of the membrane (Fig. 4A, C and S13†). The membrane headgroup region provides a favorable environment for the 1 nm particles even at high concentrations (Fig. S13†). Similar behavior can also be observed for the 3 nm size bare particle where the NPs prefer sampling the headgroup regions in the initial part of trajectories (Fig. S14†). Later on, the 3 nm NPs moved to the membrane core. Irrespective of the concentrations, the trend is adequately reproduced across systems. The NPs demonstrate characteristics considerably like partially amphiphilic small molecules, which prefer to sample the headgroup region of the lipid membranes.⁵⁹ It appears that the bare particles remained confined to the membrane (top or bottom) they initially partitioned through (Fig. S13 and S14†). The preference for this confinement was witnessed as the headgroup region for 1 nm particles and the membrane core for 3 nm particles. However, if we observe the trajectories for an additional 3 μs, we can observe a few events where the particles have reached the interface of two bilayers for B19 and B33 (Fig. S15 and S16†). Hence, transport across the multilayer skin barrier for these particles can be inferred as a rare event. As discussed earlier, the gold NPs can easily cross the skin barrier and found to be in the deeper layer. To reach into the deeper layer of the skin, the NPs must cross several layers of the lipid layers.

For the coated particles, C13 was created by initially placing two NPs in the membrane core of the top bilayer (Fig. 4B). One particle was placed in bulk water above the multilayer model (Fig. 4B). The coated NPs have distinct interactions with the membrane when compared with the bare NPs (Fig. 4B and D). It can be observed that the thiol coated particles have a preference to sample the interface of two bilayers (Fig. 4B and D), which can be observed from the trajectories of particles along the z-direction (Fig. 4D). This behavior is independent of the particle size and is found for both 1 nm and 3 nm particles (Fig. S13–S16†). However it is a strong function of the surface coating and hence was rarely observed for the bare particle as discussed previously. It is vital to sample the interface between two membranes for effective transport across membranes. Coated NPs have exhibited multiple events of sampling across the membranes, highlighting their capability to cross the multilayer skin barrier. So far, we have observed the spontaneity of gold NPs to partition into the SC membrane and how these trends get altered due to the surface characteristics of the NP. The necessity to utilize a multilayer model is noteworthy concerning its capability to study the transport across membranes, which the SLB model disregards.

The NPs cross the skin barrier via passive diffusion. Earlier studies by Yue et al.⁴⁰ and Zhang et al.⁴¹ showed that in the case of receptor mediated transport, small NPs formed clusters and translocated inside the membrane whereas bigger sized nanoparticles penetrated independently. The hydrophilic NPs’ translocation increased with an increase in their concentration but reduced for partially and totally hydrophobic NPs. In our simulations, we haven't observed any such specific translocation events; however, irrespective of membrane models and nanoparticle types, the NPs were found in the membrane interior. The coated nanoparticles formed a bigger cluster inside the membrane interior as compared to the bare NPs.

We have witnessed high headgroup fluctuations in the SLB model of the SC in the presence of NPs. Intending to determine the dependency of these fluctuations on the membrane model, we have performed a similar analysis for the DLB model. The comparison of headgroup fluctuation data between B13 and C13 for the initial 1 μs trajectory explicates an insignificant (∼0.1 nm) difference between the distributions for bare and coated NPs (Fig. 2C). Fluctuation data for the last 100 ns of the 3 μs trajectory also exhibit differences of the same order between bare and coated NPs (Fig. 2C). The fluctuation data for all the DLB demonstrate a similar trend (Fig. S17†). The lipid density (Fig. S18†) profile also confirms that even at higher concentrations of the NPs, the bilayer holds its structure and undulations were fewer as compared to the SLB. A similar conclusion can also be made based on the order parameter data for the DLB (Fig. S19†), which shows less disorder in the DLB than the SLB. Though to a lesser extent, the order parameter calculations show a higher disorder in coated particles than the bare. With these findings, we can conclude that the DLB model does not manifest high fluctuations like what was witnessed for the SLB. A multilayer or a DLB model should be employed to model the skin barrier.

Lateral density maps show a correlation between the NP density and the CHOL density for the DLB (Fig. S20–S22†), like the SLB (Fig. S4–S7†). However, the correlation is not as assertive as it was for the SLB. Dissenting to the SLB model's findings, the DLB model's density maps indicate an opposing trend. The lateral density for CER and FFA appears to be lower at the sites corresponding to the NPs. Interestingly this is observable even when the densities for lipid species are higher than that of the SLB.

Analogous to the SLB, we have calculated contacts in this model for all cases. However, there is a considerable decrease in the value of normalized contacts for the DLB compared to the single counterpart. The reduction in the case of CHO is more significant when compared with the other two species. The radial distribution function plots (Fig. S23† and S24†) also confirm these findings. Nevertheless, the relative trends are conserved from the SLB model. CHO molecules show preferential binding as compared to CER and FFA. This difference gets significantly reduced for coated NPs. Furthermore, opposite trends for preferential binding for CER and FFA in bare and coated NPs are reproduced as they were in SLB. The results highlight the selectivity of lipid species to bind to NPs based on the surface chemistry preferentially, and it is conserved in both models.

The affinity of NPs towards CHO will have a larger implication as far as the skin barrier function is concerned. The CHO is known to be an essential lipid for the better barrier function of the skin,^60–67 and due to increased interaction between the CHOL and GNP, its local concentration/orientation will change within the membrane. This reorganization may lead to local disturbance in the unique equimolar ratio of CER/CHOL/FFA responsible for strong skin barrier function. This suggests that NPs should be designed for both transdermal drug delivery, as well as their minimal impact on the skin lipid composition. This type of lipid-NP specificity opens many possibilities, such as designing NPs to use better drug delivery systems. For example, the thiol coated NPs had more contacts than bare NPs.

When both SLB and DLB simulations were analyzed for cluster formation of NPs, some common features were observed—the NPs were found in dimer, trimer and higher configurations. Fig. S25–S32† show the time evolution of NP–NP pair distances for each system. Based on these distances and visual inspection of the trajectory we found that the cluster formation process is stochastic in nature; no clear correlation was found between the cluster type (dimer, trimer, etc.) and the total number of GNPs present within a system; however a larger number of GNPs led to bigger sized clusters. Clusters formed only in the water, and they penetrated the bilayer. Post bilayer penetration no agglomeration was observed. The time required for penetration varied with respect to the cluster size. Bigger clusters penetrated slowly as compared to smaller clusters. As an example, dimers and trimers penetrated in the first 200 ns of the simulation, whereas 6-mer and higher took around 500 ns. The bigger clusters reoriented themselves to penetrate from the narrower side of the cluster (Fig. S33†).

4 Conclusion

NPs are being utilized for several personal care and healthcare applications. We have employed coarse-grained MD simulation to elucidate the interactions of bare and surface coated gold NPs with single-bilayer and double-bilayer models of the stratum corneum. Our study highlights the differences that arise in the permeation characteristics of the NPs due to the surface coating. The NPs exhibit a strong tendency to partition into the SC membrane, irrespective of their surface coating. Both SLB and DLB models provide a highly favorable environment for NPs. The SLB model gets significantly destabilized in the presence of coated NPs compared to the bare NPs as exhibited by the head group fluctuations, order parameter, and density distribution calculations. The coating of thiol groups is found to overcome the incapability of bare NPs to translocate across the membranes in the DLB. The fluctuations in the SLB seem unrealistic because they are not present in the DLB, which is a more accurate replica of the skin barrier. Our contact analyses highlight the strong tendency of lipid components of the SC layer to interact with the NPs, although CHO interacted more as compared to FFA and CER. Also, the interactions of FFA and CER with the NPs are significantly different across the bare and coated NPs.

These findings would certainly play a prominent role in targeted drug delivery across the skin barrier. The NP insertion into the SC membrane is found to be a vital function of the size and surface chemistry of the NPs. Also, targeted delivery can be employed based on the preferential interactions of the NPs with the lipid species. Our study emphasizes significant characteristics of NP–SC interactions, which highlight the necessity of effective NP synthesis for transport through the skin barriers. It also highlights the necessity of employing more realistic multi-layer models for studying the barrier offered by the SC.

Author contributions

R. G. and B. R. conceptualized the idea. Y. B., P. S. and R. G. contributed equally. All authors contributed to discussing the results and writing and editing the manuscript.

Conflicts of interest

The authors declare no conflicts of interest.

Abbreviations

NPs	Nanoparticles
SC	Stratum corneum
MD	Molecular dynamics
CHO	Cholesterol
FFA	Free fatty acids
CER	Ceramides
SLB	Single lipid bilayer
DLB	Double lipid bilayer
DPD	Dissipative particle dynamics
CG	Coarse-grained

Acknowledgements

This research was funded by Tata Consultancy Services (TCS), CTO organization. The authors would like to thank Dr Gautam Shroff, Head, TCS Research and Mr K. Ananth Krishnan, CTO and EVP, Tata Consultancy Services for their constant encouragement and support during this project.

References

1. J. M. Jungersted, L. I. Hellgren, G. B. Jemec and T. Agner, Lipids and skin barrier function – a clinical perspective, Contact Dermatitis, 2008, 58(5), 255–262.
2. P. M. Elias and J. S. Wakefield, Skin barrier function, Nutrition for Healthy Skin, 2010, pp. 35–48.
3. M. E. Lane, P. Santos, A. C. Watkinson and J. Hadgraft, Passive Skin Permeation Enhancement, Transdermal and Topical Drug Delivery: Principles and Practice, 2012.
4. G. K. Menon, New insights into skin structure: scratching the surface, Adv. Drug Delivery Rev., 2002, 54, S3–S17.
5. M. L. Clausen and T. Agner, Antimicrobial Peptides, Infections and the Skin Barrier, Curr. Probl. Dermatol., 2016, 49, 38–46.
6. J. Van Smeden and J. A. Bouwstra, Stratum Corneum Lipids: Their Role for the Skin Barrier Function in Healthy Subjects and Atopic Dermatitis Patients, Curr. Probl. Dermatol., 2016, 49, 8–26.
7. G. Cevc and U. Vierl, Nanotechnology and the transdermal route. A state of the art review and critical appraisal, J. Controlled Release, 2010, 141(3), 277–299.
8. H. Trommer and R. H. Neubert, Overcoming the stratum corneum: the modulation of skin penetration. A review, Skin Pharmacol. Physiol., 2006, 19(2), 106–121.
9. M. B. Brown, G. P. Martin, S. A. Jones and F. K. Akomeah, Dermal and transdermal drug delivery systems: current and future prospects, Drug Delivery, 2006, 13(3), 175–187.
10. A. Z. Alkilani, M. T. C. McCrudden and R. F. Donnelly, Transdermal drug delivery: innovative pharmaceutical developments based on disruption of the barrier properties of the stratum corneum, Pharmaceutics, 2015, 7(4), 438–470.
11. M. B. Brown, M. J. Traynor, G. P. Martin and F. K. Akomeah, Transdermal drug delivery systems: Skin perturbation devices, Methods Mol. Biol., 2008, 437, 119–139.
12. J. Dermol-Černe, E. Pirc and D. Miklavčič, Mechanistic view of skin electroporation–models and dosimetry for successful applications: an expert review, Expert Opin. Drug Delivery, 2020, 17(5), 689–704.
13. M. Roustit, S. Blaise and J. L. Cracowski, Trials and tribulations of skin iontophoresis in therapeutics, Br. J. Clin. Pharmacol., 2014, 77(1), 63–71.
14. J. Serup, J. Keiding, A. Fullerton, M. Gniadecka and R. Gniadecki, High-frequency Ultrasound Examination of Skin: Introduction and Guide, 2006.
15. Y. C. Kim, J. H. Park and M. R. Prausnitz, Microneedles for drug and vaccine delivery, Adv. Drug Delivery Rev., 2012, 64(14), 1547–1568.
16. P. Karande, A. Jain and S. Mitragotri, Discovery of transdermal penetration enhancers by high-throughput screening, Nat. Biotechnol., 2004, 22(2), 192–197.
17. S. Kumar, M. Zakrewsky, M. Chen, S. Menegatti, J. A. Muraski and S. Mitragotri, Peptides as skin penetration enhancers: mechanisms of action, J. Controlled Release, 2015, 199, 168–178.
18. P. Ghasemiyeh and S. Mohammadi-Samani, Potential of nanoparticles as permeation enhancers and targeted delivery options for skin: advantages and disadvantages, Drug Des., Dev. Ther., 2020, 14, 3271–3289.
19. P. Karande, A. Jain and S. Mitragotri, Insights into synergistic interactions in binary mixtures of chemical permeation enhancers for transdermal drug delivery, J. Controlled Release, 2006, 115(1), 85–93.
20. P. Karande, A. Jain, K. Ergun, V. Kispersky and S. Mitragotri, Design principles of chemical penetration enhancers for transdermal drug delivery, Proc. Natl. Acad. Sci. U. S. A., 2005, 102(13), 4688–4693.
21. R. Gupta, B. S. Dwadasi, B. Rai and S. Mitragotri, Effect of chemical permeation enhancers on skin permeability: in silico screening using molecular dynamics simulations, Sci. Rep., 2019, 9(1), 1–11.
22. N. Dragicevic, J. P. Atkinson and H. I. Maibach, "Chemical penetration enhancers: classification and mode of action, Percutaneous Penetration Enhancers Chemical Methods in Penetration Enhancement, 2015, pp. 11–27.
23. B. Kumar, S. K. Jain and S. K. Prajapati, Effect of penetration enhancer DMSO on in vitro skin permeation of acyclovir transdermal microemulsion formulation, Int. J. Drug Delivery, 2011, 3(1), 83–94.
24. R. Gupta, Y. Badhe, B. Rai and S. Mitragotri, Molecular mechanism of the skin permeation enhancing effect of ethanol: a molecular dynamics study, RSC Adv., 2020, 10(21), 12234–12248.
25. R. Notman, W. K. den Otter, M. G. Noro, W. J. Briels and J. Anwar, The permeability enhancing mechanism of DMSO in ceramide bilayers simulated by molecular dynamics, Biophys. J., 2007, 93(6), 2056–2068.
26. P. Couvreur, Nanoparticles in drug delivery: past, present and future, Adv. Drug Delivery Rev., 2013, 65(1), 21–23.
27. G. M. El Maghraby, B. W. Barry and A. C. Williams, Liposomes and skin: from drug delivery to model membranes, Eur. J. Pharm. Sci., 2008, 34(4–5), 203–222.
28. J. J. Escobar-Chavez, R. Diaz-Torres, I. M. Rodriguez-Cruz, C. L. Domínguez-Delgado, R. S. Morales, E. Ángeles-Anguiano and L. M. Melgoza-Contreras, Nanocarriers for transdermal drug delivery, Res. Rep. Transdermal Drug Delivery, 2012, 3–17.
29. Y. Huang, F. Yu, Y. S. Park, J. Wang, M. C. Shin, H. S. Chung and V. C. Yang, Co-administration of protein drugs with gold nanoparticles to enable percutaneous delivery, Biomaterials, 2010, 31(34), 9086–9091.
30. H. I. Labouta, T. Kraus, L. K. El-Khordagui and M. Schneider, Combined multiphoton imaging-pixel analysis for semiquantitation of skin penetration of gold nanoparticles, Int. J. Pharm., 2011, 413(1–2), 279–282.
31. H. I. Labouta, L. K. El-Khordagui, T. Kraus and M. Schneider, Mechanism and determinants of nanoparticle penetration through human skin, Nanoscale, 2011, 3(12), 4989–4999.
32. R. Fernandes, N. R. Smyth, O. L. Muskens, S. Nitti, A. Heuer-Jungemann, M. R. Ardern-Jones and A. G. Kanaras, Interactions of skin with gold nanoparticles of different surface charge, shape, and functionality, Small, 2015, 11(6), 713–721.
33. D. C. Liu, A. P. Raphael, D. Sundh, J. E. Grice, H. Peter Soyer, M. S. Roberts and T. W. Prow, The human stratum corneum prevents small gold nanoparticle penetration and their potential toxic metabolic consequences, J. Nanomater., 2012, 2012, 1–8.
34. D. Bedrov, G. D. Smith, H. Davande and W. Li, Passive transport of C60 fullerenes through a lipid membrane: a molecular dynamics simulation study, J. Phys. Chem. B, 2008, 112(7), 2078–2084.
35. X. Lin, Y. Li and N. Gu, Nanoparticle's size effect on its translocation across a lipid bilayer: a molecular dynamics simulation, J. Comput. Theor. Nanosci., 2010, 7(1), 269–276.
36. R. Gupta and B. Rai, Effect of size and surface charge of gold nanoparticles on their skin permeability: a molecular dynamics study, Sci. Rep., 2017, 7, 45292.
37. R. Gupta and B. Rai, Penetration of gold nanoparticles through human skin: unraveling its mechanisms at the molecular scale, J. Phys. Chem. B, 2016, 120(29), 7133–7142.
38. R. Gupta, Y. Badhe, S. Mitragotri and B. Rai, Permeation of nanoparticles across the intestinal lipid membrane: dependence on shape and surface chemistry studied through molecular simulations, Nanoscale, 2020, 12(11), 6318–6333.
39. R. Gupta and B. Rai, In silico design of nanoparticles for transdermal drug delivery application, Nanoscale, 2018, 10(10), 4940–4951.
40. T. Yue and X. Zhang, Cooperative effect in receptor-mediated endocytosis of multiple nanoparticles, ACS Nano, 2012, 6(4), 3196–3205.
41. H. Zhang, Q. Ji, C. Huang, S. Zhang, B. Yuan, K. Yang and Y.-q. Ma, Cooperative transmembrane penetration of nanoparticles, Sci. Rep., 2015, 5, 10525.
42. S. J. Marrink, A. H. De Vries and A. E. Mark, Coarse grained model for semiquantitative lipid simulations, J. Phys. Chem. B, 2004, 108(2), 750–760.
43. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. De Vries, The MARTINI force field: coarse grained model for biomolecular simulations, J. Phys. Chem. B, 2007, 111(27), 7812–7824.
44. T. A. Wassenaar, H. I. Ingólfsson, R. A. Böckmann, D. P. Tieleman and S. J. Marrink, Computational lipidomics with insane: a versatile tool for generating custom membranes for molecular simulations, J. Chem. Theory Comput., 2015, 11(5), 2144–2155.
45. J. Lin, H. Zhang, Z. Chen and Y. Zheng, Penetration of lipid membranes by gold nanoparticles: Insights into cellular uptake, cytotoxicity, and their relationship, ACS Nano, 2010, 4(9), 5421–5429.
46. J. Q. Lin, H. W. Zhang, Z. Chen, Y. G. Zheng, Z. Q. Zhang and H. F. Ye, Simulation study of aggregations of monolayer-protected gold nanoparticles in solvents, J. Phys. Chem. C, 2011, 115(39), 18991–18998.
47. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1, 19–25.
48. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, GROMACS: fast, flexible, and free, J. Comput. Chem., 2005, 26(16), 1701–1718.
49. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson and D. van der Spoel, GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit, Bioinformatics, 2013, 29(7), 845–854.
50. B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation, J. Chem. Theory Comput., 2008, 4(3), 435–447.
51. W. Humphrey, A. Dalke and K. Schulten, VMD: visual molecular dynamics, J. Mol. Graphics, 1996, 14(1), 33–38.
52. J. Wong-Ekkabut, S. Baoukina, W. Triampo, I.-M. Tang, D. P. Tieleman and L. Monticelli, Computer simulation study of fullerene translocation through lipid membranes, Nat. Nanotechnol., 2008, 3(6), 363.
53. R. Chang and J. Lee, Dynamics of C60 molecules in biological membranes: computer simulation studies, Bull. Korean Chem. Soc., 2010, 31(11), 3195–3200.
54. A. Hadad, D. L. Azevedo, E. W. Caetano, V. N. Freire, G. L. Mendonça, P. L. Neto, E. L. Albuquerque, R. Margis and C. Gottfried, Two-level adsorption of ibuprofen on C60 fullerene for transdermal delivery: classical molecular dynamics and density functional theory computations, J. Phys. Chem. C, 2011, 115(50), 24501–24511.
55. S. Zhang, Y. Mu, J. Z. Zhang and W. Xu, Effect of self-assembly of fullerene nano-particles on lipid membrane, PloS One, 2013, 8(10), 1–11.
56. Y. Guo, M. Werner, J. B. Fleury and V. A. Baulin, Unexpected Cholesterol-Induced Destabilization of Lipid Membranes near Transmembrane Carbon Nanotubes, Phys. Rev. Lett., 2020, 124(3), 038001.
57. P. Raczyński, K. Górny, J. Samios and Z. Gburski, Interaction between silicon-carbide nanotube and cholesterol domain. a molecular dynamics simulation study, J. Phys. Chem. C, 2014, 118(51), 30115–30119.
58. P. Raczyński, V. Raczyńska, K. Górny and Z. Gburski, Properties of ultrathin cholesterol and phospholipid layers surrounding silicon-carbide nanotube: MD simulations, Arch. Biochem. Biophys., 2015, 580, 22–30.
59. P. Sharma, S. Parthasarathi, N. Patil, M. Waskar, J. S. Raut, M. Puranik, K. G. Ayappa and J. K. Basu, Assessing Barriers for Antimicrobial Penetration in Complex Asymmetric Bacterial Membranes: A Case Study with Thymol, Langmuir, 2020, 36(30), 8800–8814.
60. M. Sochorová, P. Audrlická, M. Červená, A. Kováčik, M. Kopečná, L. Opálka, P. Pullmannová and K. Vávrová, Permeability and microstructure of cholesterol-depleted skin lipid membranes and human stratum corneum, J. Colloid Interface Sci., 2019, 535, 227–238.
61. E. H. Mojumdar, G. S. Gooris and J. A. Bouwstra, Phase behavior of skin lipid mixtures: the effect of cholesterol on lipid organization, Soft Matter, 2015, 11(21), 4326–4336.
62. K.-M. Joo, J.-H. Hwang, S. Bae, D.-H. Nahm, H.-S. Park, Y.-M. Ye and K.-M. Lim, Relationship of ceramide–, and free fatty acid–cholesterol ratios in the stratum corneum with skin barrier function of normal, atopic dermatitis lesional and non-lesional skins, J. Dermatol. Sci., 2015, 77(1), 71–74.
63. H. Mizushima, J. Fukasawa and T. Suzuki, Phase behavior of artificial stratum corneum lipids containing a synthetic pseudo-ceramide: a study of the function of cholesterol, J. Lipid Res., 1996, 37(2), 361–367.
64. M. Eeman, G. Francius, Y. F. Dufrêne, K. Nott, M. Paquot and M. Deleu, Effect of cholesterol and fatty acids on the molecular interactions of fengycin with stratum corneum mimicking lipid monolayers, Langmuir, 2009, 25(5), 3029–3039.
65. J. Zbytovská, M. A. Kiselev, S. S. Funari, V. M. Garamus, S. Wartewig, K. Palát and R. Neubert, Influence of cholesterol on the structure of stratum corneum lipid model membrane, Colloids Surf., A, 2008, 328(1–3), 90–99.
66. K. R. Feingold, M. Q. Man, G. K. Menon, S. S. Cho, B. E. Brown and P. M. Elias, Cholesterol synthesis is required for cutaneous barrier function in mice, J. Clin. Invest., 1990, 86(5), 1738–1745.
67. P. W. Wertz, Roles of lipids in the permeability barriers of skin and oral mucosa, Int. J. Mol. Sci., 2021, 22(10), 1–15.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2na00241h

‡ Contributed equally.

---