A Martini coarse-grained model of the calcein fluorescent dye

Calcein leakage assays are a standard experimental set-up for probing the extent of damage induced by external agents on synthetic lipid vesicles. The fluorescence signal associated with calcein release from liposomes is the signature of vesicle disruption, transient pore formation or vesicle fusion. This type of assay is widely used to test the membrane disruptive effect of biological macromolecules, such as proteins, antimicrobial peptides and RNA and is also used on synthetic nanoparticles with a polymer, metal or oxide core. Little is known about the effect that calcein and other fluorescent dyes may have on the properties of lipid bilayers, potentially altering their structure and permeability. Here we develop a coarse-grained model of calcein that is compatible with the Martini force field for lipids. We validate the model by comparing its dimerization free energy, aggregation behavior at different concentrations and interaction with a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane to those obtained at atomistic resolution. Our coarse-grained description of calcein makes it suitable for the simulation of large calcein-filled liposomes and of their interactions with external agents, allowing for a direct comparison between simulations and experimental liposome leakage assays.


Introduction
Calcein is a slightly water soluble, green fluorescent dye that is widely used to study cell viability and as indicator of lipid vesicle leakage by fluorescence microscopy imaging [1][2][3].
Calcein leakage assays are based on calcein self-quenching: at a concentration above 70 mM the fluorescence of the dye is quenched. In a typical leakage experiment, calcein is trapped into lipid vesicles or liposomes at a concentration above the self-quenching threshold. As calcein has a small permeability coefficient across a lipid bilayer [4], as long as the vesicles are intact the only fluorescence intensity that can be recorded is due to small leakage fluctuations. Varying the external conditions, such as adding an external agent to the solution, may lead to a change of the fluorescence intensity induced by calcein molecules that have leaked out from the vesicles to the surrounding environment, where calcein is dissolved at a smaller concentration. This is a signature of vesicle disruption, pore formation or vesicle fusion.
This type of assay is widely used to test the membrane disruptive effect of proteins [1], antimicrobial peptides [5], poly-cations [6], single strands of RNA [2] and also nanoparticles [7][8][9] (NPs). Inorganic NPs, such as metal [10,11] or oxide [12,13] NPs, are nowadays considered as multivalent agents for biomedical applications, ranging from diagnostics to therapy. Many of these applications would require some degree of control on NP-membrane interactions: once the NPs have reached their cell target, non-disruptive membrane translocation or transient membrane damage can be effectively exploited to deliver drugs to the cytoplasm, while the disruption of significant portions of the plasma membrane can be lethal to the cell.
In this perspective, liposome leakage assays, possibly performed on synthetic vesicles, are fundamental to rationalize which physico-chemical factors drive the NP towards a gentle or disruptive type of interaction with the lipid bilayer. The role played by the NP charge, for example, has been investigated using calcein leakage assays by Goodman [8]  Similarly, Liu [14] et al. has shown the transient poration due to the interaction of zwitterionic liposomes with citrate-capped Au NPs. Calcein leakage assays have been used to test the effect of varying the NP core material [15], size, shape, and surface functionalization [16] as well.
The computational approach, and the one based on molecular dynamics (MD) in particular, can certainly be a useful complement to the experiments performed on synthetic vesicles.
Simulations offer the unique opportunity to track NP-membrane interactions with atomistic or molecular details. A fruitful comparison between experimental and in silico results, though, relies also on the possibility to simulate systems whose composition is as close as possible to that of the experimental samples. Fluorescent probes, then, come into play.
Little is known about the effect that the many different fluorescent dyes used in the experiments, sometimes at very large concentration, may have on the properties of lipid bilayers or on the interactions between them and the NPs. Calcein is not an exception: at physiological pH the molecule is expected to bring four negatively charged COOgroups, which may lead to a favorable interaction with the polar lipid head region. At large concentrations (> 100 mM) this interaction may have effects on the structural and mechanical properties of the bilayer, in turn affecting bilayer permeability and its interaction between external agents or NPs. For these reasons, it is of interest to study the calceinmembrane interactions using a computational approach.
Here we present the development of a coarse-grained (CG) model of calcein compatible with the popular Martini force field for lipids [17][18][19]. The model resolution makes it suitable to the simulation of calcein-filled liposomes and of their interactions with external agents. We base our model on the structural properties of calcein, as obtained from atomistic simulations, and validate it by comparing the aggregation properties of the CG calcein in water to those obtained at atomistic level.
As calcein is charged at physiological pH, it is interesting to evaluate the performance of different versions of the Martini force field, which treat differently the electrostatic interactions. We compare three versions of Martini: the first is the standard Martini forcefield [17] (STD). In STD Martini, electrostatic interactions are described as short-range interactions and the water polarizability is taken into account implicitly with a uniform dielectric constant er = 15. The long-range contribution to the electrostatic interactions can be added to the STD model via the particle-mesh Ewald (PME) method, an approach that has been shown [20] to improve the performance of the model when dealing with highly charged molecules. Here, we will always include the PME energy term in our STD simulations. The second model we use includes explicitly water polarizability [18] (PW). The model uses PME and a uniform dielectric constant er = 2.5. We also test the polarizable water model (refPW) developed by Michalowsky [19] et al., in which water dielectric properties are refined with respect to the original PW model.
In the Methods section we provide all the details about the force field parameters, MD run set-up and system composition; in the Results section we present the development of the CG Martini model of calcein, and its validation in terms of dimerization free energy, aggregation behavior at different concentrations, and interaction with a 1-palmitoyl-2-oleoylsn-glycero-3-phosphocholine (POPC) membrane. Eventually, the different CG model performances and perspective used are discussed.

Methods
Both atomistic and CG simulations, as well as the analysis, are performed with the Gromacs 2016 molecular dynamics software [21].

Atomistic simulations
Atomistic simulations are performed with the 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] (tT = 2 ps); the pressure is kept constant at 1 bar by the Berendsen barostat (tp = 1 ps) for the equilibration runs and by the Parrinello-Rahman algorithm [26] (tp = 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 5x5x5 nm 3  Simulations are initialized with the calcein molecule placed on top of the membrane, in the water phase.

Coarse-grained simulations
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] (tT = 1 ps), Berendsen

Umbrella sampling simulations
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 production runs are 500 ns long. The total simulated time is about 1.5 µs (atomistic) and 13.7 µs (CG).
For what concerns the umbrella sampling simulations of the desorption process of one calcein molecule from the POPC membrane we use, as order parameter, the projection along 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 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).

Contacts and aggregates analysis
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 [30] et al. 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.

Atomistic model
The parameterization of the atomistic model of calcein is based on the OPLS-UA force-field developed by Jorgensen [22]   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].
We added improper dihedrals to impose the planarity of the four ethanoate groups (-CH2-COO -), that of the three rings of the central region, that of the g-butyrolactone and that of the benzene rings. In the case of the two sp 3 N atoms, we added an improper dihedral in order to keep the tetrahedral geometry. The parameterization of the benzene ring is the one proposed by De Jong [39] et al.
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 kJmol -1 rad -2 ) 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.

Validation
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.
Dimerization process. We calculate the potential of mean force (PMF) of the dimerization of two calcein molecules in 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 the details about the simulation set-up see Section 2. Each profile is also shifted by the factor 2kBT ln(d) in order to take into account the Jacobian of the transformation from Cartesian to spherical coordinates [42].
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.     For what concerns 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 1phthalanone face the interior of the micelle. On the contrary, the 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 better reproduces the atomistic aggregation behavior, characterized by the formation of large aggregates with possible detachments on time scales of microseconds;  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 Figure 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 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 interacts with the phosphate groups only, and can be found at a distance of 1.2 nm (STD) and 0.7 nm (atomistic). The simulation with the PW model, instead, 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 also, with small probability compared to the STD and atomistic models, it predicts a bound metastable state at a distance of 0.25 nm from the phosphate group (see Figure S8).
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

Discussion
We can now sum up and discuss the similarities and the differences between the three CG models of calcein we have tested. All models predict that the dimerization process is thermodynamically favored (Figure 3 interactions that has been previously 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 p-p 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 confirmed also by the unbiased simulations of calcein aggregation. The number of ring-ring contacts is overestimated by all CG models, and especially by the STD model. In terms of aggregation, the agreement between the PW models and the atomistic model, though, 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. At atomistic resolution, on the contrary, 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.
The Nda-Nda interactions, though, are stronger than the Na-Na interactions (4.5 kJ/mol vs. 4 kJ/mol), and thus provide a better aggregation behavior. We have tried also to increase further 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 especially in terms of aggregation at high concentration, the behavior of the STD (with long-range electrostatics) model seems in qualitative agreement with the atomistic one.
Experimentally, fluorescence-lifetime measurements [41] suggest that calcein's selfquenching 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; PW models, instead, 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). The PW model, on the contrary, 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 of 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 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, especially when looking at calcein in 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 to include the leaking dye in the simulations may help to clarify the role played by the dye itself during membrane poration.
Calcein is used also as model metabolite to study the permeability of gap junction channels.
Zonta [47] et al., for example, used calcein to probe the permeability of the human connexin26 channel protein to negatively charged molecules bearing different total charge. 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 assess its role during translocation [48], accessing longer time scales than with a pure atomistic approach.

Conclusions
In this paper we have developed and tested a coarse-grained model of calcein, a fluorescent dye routinely used in liposome leakage assays. The model is compatible with the Martini force field, and as such is suitable to the simulation, via molecular dynamics, of the interaction between synthetic nanoparticles and large lipid bilayers (e.g. calcein-containing liposomes with diameters of tens of nm). The availability of a coarse-grained model of calcein is a step forward in the direction of a better integration of experimental and computational data, as the model allows to perform simulations in which the effect of the dye during the membrane poration process is taken into account explicitly and with molecular resolution.