A Martini coarse-grained model of the calcein fluorescent dye

Atomistic simulations are performed with the OPLS united-atom (OPLS-UA) force field [22], the TIP3P water model [23] and the Berger parameters for lipids [24]. The time step is set to 2 fs. All simulations are performed in the NpT ensemble: the temperature is kept constant at 310 K by the velocity rescale thermostat [25] (τ_T  =  2 ps); the pressure is kept constant at 1 bar by the Berendsen barostat (τ_p  =  1 ps) for the equilibration runs and by the Parrinello–Rahman algorithm [26] (τ_p  =  1 ps) for the production runs. The cut-off of the Van der Waals interactions is set to 1 nm, while the electrostatic interaction is treated via the PME method with a grid spacing of 0.12 nm.

The simulation of a single calcein in water is made in a cubic box of 5  ×  5  ×  5 nm³ filled with water and four Na⁺ counterions. The unbiased simulations with calcein at different concentrations are made in a cubic box of 10  ×  10  ×  10 nm³ containing water, salt at physiological concentration and Na⁺ counterions. At 30 mM calcein concentration, 20 calcein molecules are solvated by about 32 000 water molecules; at 150 mM calcein concentration, 100 calcein molecules are solvated by about 30 000 water molecules. The simulations for the dimerization process are made with two calcein molecules solvated by about 33 000 water molecules and counterions in a dodecahedron box (radius of the inscribed sphere about 5.5 nm). To study the interaction between calcein and a model zwitterionic lipid membrane we set-up a cubic box of 6  ×  6  ×  10 nm³ with a symmetric membrane composed of 114 POPC lipids, 7100 water molecules, salt at physiological concentration and counterions. Simulations are initialized with the calcein molecule placed on top of the membrane, in the water phase.

The CG simulations are performed with the Martini force field [17]. The time step is set to 20 fs. All simulations are performed in the NpT ensemble with temperature and pressure set to 310 K and 1 bar, respectively. The velocity rescale thermostat [25] (τ_T  =  1 ps), Berendsen (τ_p  =  4 ps) and Parrinello–Rahman [26] (τ_p  =  12 ps) barostat are used. For the simulations with the STD model, the cut-off of the Van der Waals interactions is set to 1.1 nm and the dielectric constant is ε_r  =  15. The simulations performed with the PW model have a Van der Waals cut-off set to 1.2 nm and a dielectric constant ε_r  =  2.5. For the refPW model the Van de Waals cut-off is 1.1 nm. In all cases, the Verlet scheme for the neighbor list update and the PME method for the electrostatic interactions, with a grid spacing of 0.12 nm, are used.

The simulation boxes are set-up in the same way as in the atomistic case with the same box dimensions and number of calcein molecules. A factor of 4 to 1 is considered to rescale the number of water molecules from atomistic to CG. Counterions are always present and ions at physiological concentration are present as in the atomistic case.

The analyses of the umbrella sampling [27] simulations are performed with the WHAM [28] algorithm using the 'gmx wham' tool [29] of the Gromacs suite. The order parameter for the umbrella sampling of the dimerization process of two calcein molecules is the distance between their centers of mass (COMs). We sample, both at atomistic and CG level, 25 windows: from a distance of 0.4 nm to 5 nm. In the atomistic case, each window is equilibrated for 10 ns followed by a production run of 50 ns. In the CG case, the equilibration runs are 50 ns long while the production runs are 500 ns long. The total simulated time is about 1.5 µs (atomistic) and 13.7 µs (CG).

With regards to the umbrella sampling simulations of the desorption process of one calcein molecule from the POPC membrane we use, as an order parameter, the projection along the z-axis of the distance between the COM of the membrane and the COM of the calcein molecule. In this case 17 windows are sampled: from 1.6 nm to 4 nm. Each window is initialized with a frame extracted from the corresponding unbiased simulation. At about 1.6 nm the calcein molecule is partially embedded in the membrane head region while, at a distance larger than 5 nm, it is in the water phase. Each window of equilibration and production runs lasts 15 ns  +  75 ns (atomistic) and 50 ns  +  800 ns (CG). The total simulated time is about 1.5 µs (atomistic) and 14.5 µs (CG).

The contacts analysis is performed with the 'mindist' Gromacs tool. Two atoms or beads are considered in contact if their distance is lower than a threshold of 0.8 nm. The analysis of the aggregates is performed with the 'g_aggregate' tool, a Gromacs compatible tool developed by Barnoud et al [30]. Two calcein molecules are assigned to the same cluster if the minimum distance between their atoms or beads is lower than a cut-off value which we set to 1.2 nm.

The parameterization of the atomistic model of calcein is based on the OPLS-UA force field developed by Jorgensen et al [22]. The full details of the topology are available in our open online repository [31].

In the left panel of figure 1 the chemical structure of calcein and its atomistic UA representation are shown. We have a central hydrophobic region that is highly rigid and planar and the 1-phthalanone perpendicular to the central region. The atom type assignment, the bonded and non-bonded interaction parameters, as well as the partial charge assignment, are derived from the OPLS-UA parameters for hydrocarbons and proteins [22, 32, 33] and from the AMBER parameters for nucleic acid and proteins [34, 35].

Zoom In Zoom Out Reset image size

We added improper dihedrals to impose the planarity of the four ethanoate groups (–CH₂–COO⁻), that of the three rings of the central region, that of the γ-butyrolactone and that of the benzene rings. In the case of the two sp³ N atoms, we added an improper dihedral to keep the tetrahedral geometry.

3.2.1. Mapping.

3.2.2. Bonded interactions.

The parameterization of the CG bonded interactions is based on the reproduction of the target atomistic distributions of distances and angles. Atomistic distributions are obtained from a trajectory of a single calcein molecule in water [38]. In figures S1 and S2 (stacks.iop.org/JPhysD/51/384002/mmedia) of the supplementary material, bond and angle distributions are shown, respectively. The equilibrium distances or angles and the force constants of the harmonic (or cosine-harmonic for angles) functions are chosen so as to reproduce the peak position and the width of the target distributions. Bond constraints are used whenever the force constant needed to reproduce the atomistic distribution is larger than 12 000 kJ mol⁻¹ nm⁻². The parameterization of the benzene ring is that proposed by De Jong et al [39].

In order to keep the central region of the molecule planar, we used two improper dihedrals (between beads 4–6–5–8 and 7–8–5–6, see figure 1) with equilibrium angles and force constants chosen to fit the atomistic distributions, as shown in figure S3. To keep bead 9 in the same plane of the central region we added an improper dihedral between beads 5–6–8–9 and the corresponding angle distribution is shown in the left panel of figure 2. The position of the 1-phthalanone is kept perpendicular to the central region by an improper dihedral (beads 9–12–10–11, equilibrium angle of 180° and force constant of 80 kJ mol⁻¹ rad⁻²) and by an angle interaction between beads 5–9–10. At CG level, beads 9–12–10–11 are coplanar, as shown in the right panel of figure 2.

Zoom In Zoom Out Reset image size

In calcein release experiments the calcein solution inside the liposome, at high concentration (>70 mM), is quenched by both static and dynamic quenching mechanisms [40, 41]. We thus validate our CG models based on dimerization free energies and aggregation behavior of calcein at different concentrations.

3.3.1. Dimerization process.

We calculate the potential of mean force (PMF) of the dimerization of two calcein molecules in a water solution. The PMF is obtained from umbrella sampling simulations in which the order parameter is the distance between the COM of the two molecules. For details of the simulation set-up see section 2.

A comparison of the PMF between the atomistic and CG models is shown in figure 3. As we can see, all models agree that the bound state is more stable than the dissolved state. The free energy difference between the unbound and bound state is about 17 kJ mol⁻¹ for the atomistic model. The STD and refPW model overestimates it, with a free energy difference of ~29 kJ mol⁻¹ and ~20 kJ mol⁻¹, respectively, while the PW model slightly underestimates it, with a difference of ~11 kJ mol⁻¹. The atomistic, PW and refPW models predict a small barrier for aggregation, that is not captured by the STD model. In all cases, the free energy barriers for dimerization and separation are of a few k_BT, allowing for the dynamic formation and disruption of small aggregates on experimentally relevant time scales. The width of the dimer state well is not the same in the atomistic and CG models: the latter (particularly the PW model) shows two dimer states that slightly differ in energy. In figure 4 the two states are shown. The CG dimer that corresponds to a distance between the calcein molecules of about 0.5 nm (left panel of figure 4, state A) is the same as that observed in the atomistic case. State A is more populated in the CG simulations with the STD model. This state is characterized by the flat contact between the central planar regions of calcein, with the 1-phthalanone ring pointing to the outer direction. The other dimer state, at about 0.8 nm (right panel of figure 4, state B), is the one in which the benzene rings are facing each other in a sandwich π–π stacking. This state is the most populated in the CG simulations with the PW and refPW models and is not observed in atomistic simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3.2. Aggregation behavior.

The self-quenching concentration threshold for calcein is about 70 mM. We set-up atomistic and CG simulations at two different calcein concentrations, well above and well below the threshold: 30 mM and 150 mM. We use the atomistic system as a reference and compare it to the STD and PW systems at both concentrations; for the system at 150 mM we also test the refPW. Figure 5 shows the number of aggregates as a function of time as obtained with the different models.

Zoom In Zoom Out Reset image size

At both concentrations, but with different kinetics, the atomistic simulations predict that the calcein molecules aggregate into one single cluster. The STD model agrees with the atomistic result at both concentrations. On the contrary, the CG models with water polarizability are not able to reproduce this effect in any of our simulations. By visualizing the PW or refPW trajectories, we observe a dynamic formation-destruction of micelles of calcein molecules that do not tend to form a stable cluster. The time scales of the simulations are certainly too short to exclude the possibility of disaggregation events for the atomistic and STD models: indeed, at large calcein concentration, both the atomistic and the STD CG model predict transient disaggregation events, on time scales of microseconds.

With regards to the structure of the aggregates, we also observed that the PW model has the tendency to stabilize micellar aggregates in which the benzene rings of the 1-phthalanone face the interior of the micelle. On the contrary, visual inspection of the atomistic runs suggests that most of the benzene rings are in the water phase. We calculated the number of benzene–benzene contacts, as shown in figure 6. We see that the CG models give a number of benzene–benzene contacts larger than the atomistic model by a factor of 1.5–2. We conclude that the STD model is, among the CG models tested here, the one that best reproduces the atomistic aggregation behavior, characterized by the formation of large aggregates with possible detachments on time scales of microseconds; as for the structure of the aggregates, all the CG models tested share the limitation of an overestimation of the benzene–benzene contacts.

Zoom In Zoom Out Reset image size

3.3.3. Calcein-membrane interactions.

To test the calcein-membrane interactions at atomistic and CG level, we ran a set of unbiased simulations (5.7 µs at atomistic level and 20 µs at CG level, for each model) comprising a POPC membrane (114 lipids) and a single calcein molecule, initially placed far on top of the membrane surface.

For the atomistic and STD models the calcein molecule binds to the membrane and unbinds spontaneously. Binding and unbinding events can be monitored by tracking the distance between calcein and the COM of phosphate groups, as shown in figures S6 and S8. Two metastable states can be distinguished in which the calcein molecule interacts with the membrane. In the first state, or bound state, the hydrophobic benzene ring of calcein interacts with the glycerol groups of the lipids. The correlation between the calcein-COM of membrane distance and the number of contacts between the benzene ring of calcein and the lipid glycerol groups is shown in figure S7. A snapshot of this membrane-bound configuration is shown in figure 7 for both the atomistic and STD resolutions. In CG simulations, in this configuration the calcein molecule is at a distance of about 0.1 nm from the phosphate groups; in atomistic simulations the distance is about 0.2 nm. In the second metastable state, or adsorbed state, the calcein molecule only interacts with the phosphate groups, and can be found at a distance of 1.2 nm (STD) and 0.7 nm (atomistic). On the contrary, the simulation with the PW model does not sample any bound or adsorbed state in the simulation time. The refPW model still predicts that the water phase is the most favorable state, but with small probability compared to the STD and atomistic models, it also predicts a bound metastable state at a distance of 0.25 nm from the phosphate group (see figure S8).

Zoom In Zoom Out Reset image size

To be more quantitative, we calculated the PMF of the binding process of one calcein from the water phase to the membrane-bound state, showed in figure 8. The energy profile reveals two metastable states for both the atomistic and the STD model. The refPW model predicts one metastable bound state, but its most favorable configuration is in the water phase. The PMF confirms the absence of any stable state near the membrane for the PW model. The STD energy barrier to go from the adsorbed state to the bound state is about 3.5 kJ mol⁻¹, while the adsorbed state is more stable than the bound state by ~1 kJ mol⁻¹. The barrier for the atomistic model is higher, and it sets to 5.4 kJ mol⁻¹. The energy barrier connecting the dissolved state to the bound state for the refPW model is 6.4 kJ mol⁻¹.

Zoom In Zoom Out Reset image size

We can now summarize and discuss the similarities and the differences between the three CG models of calcein we have tested. All the models predict that the dimerization process is thermodynamically favored (figure 3), in agreement with the prediction of the atomistic force field. Only the polarizable water models (PW and refPW) are able to reproduce the energy barriers of the binding and unbinding processes. The STD model overestimates the population of the dimer state over the dissolved one and does not describe the small barrier for aggregation that is foreseen by the PW, refPW and atomistic models. The behavior of the STD model, which we have employed here with the addition of long-range electrostatic interactions, is likely to be attributed to an overestimation of hydrophobic interactions that has previously been reported for other aromatic compounds [43, 44]. A more detailed inspection of the CG dimers reveals the existence of a metastable minimum in which the benzene rings of the 1-phthalanone are in π–π sandwich stacking, a configuration that is not sampled during the atomistic umbrella sampling calculations.

The overestimation of the hydrophobic interactions by the CG models is also confirmed by the unbiased simulations of calcein aggregation. The number of ring-ring contacts is overestimated by all CG models, particularly by the STD model. However, in terms of aggregation, the agreement between the PW models and the atomistic model is not as satisfying as one could expect based on the dimerization free energy profile. Both the CG models that include water polarizability (PW and refPW) tend to form a highly dynamic number of calcein micelles that do not aggregate in a stable cluster. These micelles are formed by a variable number of calcein molecules, depending on the calcein concentration, as shown in figure S5 in which the maximum size of the micelles is plotted as a function of the simulation time. A common characteristic of these micelles is that the packing of the molecules is driven by the 1-phthalanone rings that tend to stick to each other in the inner part of the micelle. On the contrary, at atomistic resolution the packing of the molecules is driven by the stacking of the central plane of calcein. During the model optimization we explored alternative mappings with the aim to improve the performance of the PW models in this respect. We have tried to have a less hydrophobic benzene with the introduction of a new Martini bead type that is equal to the SC5 bead but has a decreased self-interaction (−6%) and an increased interaction with the water beads (+6.5%). The improvement is little and does not justify the introduction of a new bead type. The final PW model uses an Nda type to describe the central C–O–C group of calcein. A straightforward choice based on standard Martini mappings would suggest to use the Na type, which has only hydrogen acceptor character. However, the Nda–Nda interactions are stronger than the Na–Na interactions (4.5 kJ mol⁻¹ versus 4 kJ mol⁻¹), and thus provide a better aggregation behavior. We have also tried to further increase the self-interaction of the beads of the central plane (beads 4, 5 and 7, up to  +11%). Again, this fine tuning did not significantly improve the aggregation behavior of the PW or refPW models. Overall, and particularly in terms of aggregation at high concentration, the behavior of the STD (with long-range electrostatics) model seems in qualitative agreement with the atomistic model.

Experimentally, fluorescence lifetime measurements [41] suggest that calcein's self-quenching method is based on both static quenching and dynamic quenching, the latter prevailing over the former. The results obtained with the atomistic and CG STD model would suggest that the static quenching deriving from calcein-calcein aggregation (over time scales longer than the fluorescence lifetime, 4 ns) is the main mechanism of fluorescence quenching; on the other hand, PW models show an aggregation behavior which is more consistent with a dynamic quenching mechanism.

Our simulations show a good agreement between the STD model predictions and the atomistic ones in terms of interaction with the lipid membrane. Both the STD and atomistic models indicate the existence of two metastable states in which the calcein is in close contact with the membrane (the bound state is shown in figure 7). On the contrary, the PW model shows no interaction at all between the negatively charged calcein and the membrane surface (figure 8). The refPW model slightly improves the situation, but it does not predict stable binding either (figure 8). This latter behavior is at odds with the general experimental observation of stable interactions between phosphatidylcholine membranes and anionic macromolecules [45] and nanoparticles [9, 46]. A possible explanation for this behavior could be the much larger hydration of the lipid head groups that is provided by the PW model. We calculated the density of water as a function of the z component of the distance from the COM of a POPC membrane, in the absence of calcein, with the OPLS-UA, STD, PW and refPW models. Table S1 contains the water density recorded in the lipid phosphate plane and shows that the STD underestimates the atomistic water density by 4%, while PW and refPW overestimate it by 11% and 14%, respectively. Head group hydration could screen the favorable interaction of the negatively charged groups of calcein with the membrane choline groups. These observations lead us to the conclusion that the use of the STD model of calcein could be preferred over the use of the PW model, particularly when looking at calcein in the presence of a lipid membrane.

We anticipate that our CG model of calcein will be useful to the interpretation of vesicle and liposome leakage assays. Different external agents, such as cell penetrating peptides, polymer and inorganic NPs, can damage the lipid bilayer via different mechanisms and to a different extent. In order to achieve a molecular interpretation of the experimental data, the possibility for inclusion of the leaking dye in the simulations may help to clarify the role played by the dye itself during membrane poration.

Calcein is also used as a model metabolite to study the permeability of gap junction channels. Zonta et al [47], for example, used calcein to probe the permeability of the human connexin 26 channel protein to negatively charged molecules bearing different total charges. They found that while the neutral form of calcein is able to cross the channel and the transition rate is comparable to that obtained experimentally, the channel is impermeable to the fully deprotonated (−4e⁻) calcein. With our CG model it would be relatively easy to change the protonation state of the dye and to assess its role during translocation [48], accessing longer time scales than with a pure atomistic approach.