Excimer formation of pyrene labeled lipids interpreted by means of coarse-grained molecular dynamics simulations

The excimer formation dynamics of pyrene-labeled molecules in lipid bilayers depends on molecular motion over distances of the order of 1-2 nm. From the concentration dependence of the excimer photoemission curve, it is possible to derive a value for the lipid self-diffusion coefficient. This technique has been intensively used in the past twenty years, leading to rather large numerical values for self-diffusion compared with other approaches based on fluorescent probes tracking. In most cases, the interpretation of the experimental data rely on models for diffusion limited 2d reaction rates, or comparison with 2d lattice random walks. Our approach uses realistic molecular dynamics trajectories to reinterpret these experiments. Based on a well established coarse-grained model for lipid MD simulations (Martini), we show how to relate simulation results to experimental data on excimer formation. Our procedure is quite general and is applicable to all diffusion-limited kinetic processes. Key to our approach is the determination of the acceleration factor of lipid coarse-grained numerical models compared to reality. We find a significant reduction of the diffusion coefficient values, in particular when interleaflet association is taken into account. Our work does not point to deviation from a diffusion-limited mechanism but indicates that the excimer formation across bilayer leaflets could be hindered.


I. INTRODUCTION
Most lipid self-diffusion determination methods rely on long range diffusion (Fluorescence Recovery After Photobleaching FRAP, Fluorescence Correlation Spectroscopy FCS, Pulse Field Gradient NMR). On the other edge, inelastic neutron scattering probes the lipid dynamics on very short time and length scales. The pyrene excimer formation approach, initiated long ago [1], is unique as far its characteristic time (100 ns) and length (1-2 nm) scales are concerned. It opens a window onto the nanometric scale in lipid membrane organization, currently subject to many supposed phenomena (rafts, static or dynamic nanodomains). How much is lipid diffusion dynamics on these scale regular or anomalous is a topic of great interest.
So far, numerical simulations are promising investigation tools to answer these questions. On the other hand, pyrene excimer formation experiments have been mostly interpreted in terms of 2d lattice random walks properties, which cannot describe faithfully the correlated motion of lipid molecules at this scale. There is little doubt that an analysis based on realistic molecular motion should outperform any approach based on point-like particles jumping over fixed discrete lattices.
In this article, we attempt to refine the past analysis of pyrene-lipid analogues diffusion by using coarse-grained (CG) molecular dynamics (MD) based on the celebrated Martini force-field [2]. Introducing a novel statistical approach that relates molecular dynamics trajectories to diffusion limited association kinetics among lipid chain subgroups, we obtain a numerical analogue to the excimer/monomer emission ratio as a function of the probe concentration. This allows us to obtain the experimental diffusion coefficient by matching the short time excimer formation rate to the long time diffusion displacement.
We find diffusion coefficients values significantly smaller than previously admitted. This provides an important clue as far as reducing the spread among all experimental values currently found in the literature.
Our paper is organized as follows. Section II introduces the principle of excimer formation dynamics. Section III explains how one can relate MD trajectories to experimentally relevant data. In Section IV we analyze data published by Vauhkonen et al., and Sassaroli et al. [3,4]. Results are discussed in Section V followed by Conclusion.

II. EXCIMER FORMATION DYNAMICS
Excimers are complexes formed by two identical molecules, one being in an excited state, the other in its ground state. The fluorescence emission spectrum of the isolated molecules, commonly referred as monomers, markedly differs from the emission spectrum of the excimers, or dimers. This makes it possible to measure optically the fraction of fluorescent probes ending up into an excited dimer (excimer) state.
Excimer formation results from a diffusion limited mechanism of bimolecular reaction kinetics, and the formation rate is related to the self-diffusion coefficient of the fluorescent probes [5]. The probability that an excited monomer ends up forming an excimer is the outcome of a competitive process involving bimolecular collision rate on the one hand, and spontaneous de-excitation rate (inverse fluorescence lifetime) on the other hand.
Pyrene molecules form very stable excited dimer states. They can be incorporated, as pyrenyl groups, into phospholipid analogues and inserted in lipid mem-branes. These fluorophores are characterized by a longer than usual lifetime, ranging from 100 to 200 ns. The pyrene excimer formation is therefore reporting on the dynamical moves of lipid molecules on this time scale, corresponding to spatial displacements of the order of 1 to 2 nm in the bilayer. In this respect, this dynamics probes local molecular motions at a nanometric scale. The use of pyrene derivatives for the purpose of studying lipid membrane dynamics was reviewed by Somerharju [6]. A larger class of diffusion-limited excitation or deexcitation processes in membranes is covered in a review by Melo and Martins [7].
The determination of lipid diffusion coefficients D using pyrene derivatives is quite indirect.
Fluorimetric techniques provides a "titration curve" of the monomer and excimer normalized emission intensities J M (x), J E (x) as functions of the pyrene molar fraction x and monomer fluorescence lifetime τ M . Statistical models of molecular diffusion must be used to relate the diffusion constant to the excimer formation rate. So far, approaches based on cubic and hexagonal lattice random walks, and Smoluchovsky-Naqvi theoretical models have been used to interpret the experimental data [1,3,4,[8][9][10][11]. The approaches based on lattice random walks tend to give larger diffusion coefficient values than the ones obtained, e.g., by tracking fluorescent lipid probes over larger scales, such as in FRAP or FCS experiments (see Discussion).
Current theoretical approaches are limited to ideal random walks or Wiener diffusion processes and they also involves approximations [12][13][14][15]. Anomalous Brownian displacements of the pyrene groups, membrane heterogeneous nanodomains could both invalidate the current numerical estimates of lipid diffusion coefficients. Conversely, a robust microscopic description of pyrenyl derivatives excimer formation would certainly help to improve the description of short range dynamical properties of lipid membranes.
In this work, we propose to use a well-known coarsegrained model for lipid bilayers, the Martini force field, and to analyze the excimer fluorimetric data based on realistic molecular dynamics (MD) trajectories.

A. First-passage reaction kinetics
Diffusion-limited bimolecular reactions represent a class of physical or chemical kinetics in which reagents react at first encounter with a yield equal to, or close to unity. This limit means that there is no activation energy barrier, nor orientation barrier opposing the formation of the product [16].
The reaction rate of a diffusion limited reaction depends exclusively on the diffusion dynamics of the two species and on the geometrical shape of the reaction region. The most natural treatment, initiated by Smolu-chovsky, considers two spheres following ideal Brownian trajectories and reacting as soon as their relative distance drops below a critical value: the capture radius ρ c . Naqvi derived a Smoluchowsky rate in the two dimensional case [12]. Both approaches predict a timedependent reaction rate and are restricted to ordinary (Wiener) diffusion processes. They are mean-field in nature as relative spatial correlations of the reactants are neglected [14].
The experimental results reported in [4] were obtained with a pyrenyl group attached to a 10 atoms long alkyl chain. This group is mobile with respect to its parent molecule center of mass, and its Brownian displacement is affected by the intramolecular degrees of freedom on short time-scales. Its motion cannot be represented accurately by a Wiener process, as the associated meansquared displacement is not just a linear function of time and includes intramolecular Rouse dynamics. The consequences of such deviations with respect to the simplest Brownian model description, on the collisional dynamics properties of the pyrenyl groups must be assessed.
The Martini force-field treats each phospholipid molecule as an assembly of 10 to 14 beads, according to a 4 atoms to 1 bead correspondence. This level of coarsegraining preserves to a certain extent the intramolecular conformational degrees of freedom of the original atomistic model. With access to submolecular details, one can attempt to improve the collisional dynamics description compared with lattice or continuous models reducing lipids to point-like objects. This is the rationale behind our treatment of pyrene derivative excimer formation dynamics.
On the other hand, CG models are characterized by a nonphysical accelerated time scale, so that they cannot predict transport properties without a reference element (experiment or all-atoms simulation). The reference element, in our case, is the experimental monomer/excimer fluorescence ratio.

B. Neutrality of fluorescent lipid probes substitution
Fluorescence is an extremely useful and versatile technique. Its use in the field of membrane studies includes imaging, life-time, anisotropy depolarization and resonant energy transfer studies. In most cases, the amount of necessary fluorophores is small, as fluorescence detection can be very sensitive. A common assumption of all these approaches is that the fluorescent labeled molecules, which often are made as similar to phospholipid molecules as possible, modifies only marginally the structure and the dynamics of the host membrane. For instance, the diffusion coefficient obtained in a FRAP experiment, is the one associated with the fluorescent probes. The identity between probe and normal lipid motions remains an assumption.
The pyrene excimer formation method requires rela-tively high amounts of pyrenyl derivatives, up to 10% molar ratio. It is assumed that the motion of the pyrenyl-grafted lipids is identical to the regular lipid compounds, and that the bilayer structure is not modified significantly. We therefore adopt this view and treat pyrene-labeled lipids as if they were ordinary phospholipid molecules. This is both a simplifying assumption and a limit of the present approach. We do not have at our disposal a coarse-grained model of pyrene-labeled lipid molecule. The polyaromatic pyrenyl group is clearly bulky compared with standard alkyl chains, and the effective van-der-Waals interaction parameters should be increased to reflect its higher polarizability. It is difficult to estimate these parameters at coarse-grained level, and guesswork is unlikely to provide an acceptable set of interactions. Another difficulty arises when it comes to model the interaction between a photoexcited pyrene group and one in its ground state. It is believed that attractive interactions arise from resonant energy states, leading to stronger and longer range ( r −3 with separation r) attractive forces, eventually leading to the complexed dimer state [5,17]. In a realistic model, those modified interactions should be estimated and taken into account.
Even if a CG model of pyrenyl derivatives existed, simulating a diluted mixture of probes and repeated sequences of excitation → diffusion → excimer formation would require lengthy and multiple runs, to obtain just a rough statistics of the events. Our opinion is that such an intense computational effort would only make sense if a trusted and robust CG model of pyrenyl compounds was first established, including the modified interactions of the photoexcited pyrene group.
We remain faithful to the implicit assumption that fluorescent probes behave the same as the major lipid components. Ordinary MD simulation runs are performed and trajectories acquired. The equilibrated trajectories are then analyzed a posteriori, as if a few of the lipid beads were actually pyrenyl moieties. Our strategy is therefore to rerun the trajectory (without recomputing the forces) with a relabeling of some beads a posteriori. In practice, one can reduce the determination of the pyrene excimer formation assays to regular dynamical time dependent correlation functions of the ordinary lipid molecules. The next subsection introduces such dynamical quantities.

C. Survival probabilities
A CG model of DOPC is represented on Fig. 1 and Fig. 2. We propose to assign to the pyrene group a set of three beads numbered from 12 to 14, hereafter referred as the the fluorescent group, or terminal subgroup (TS). A similar procedure defines the POPC terminal subgroup (not shown). We assume that the excimer association takes place as soon as the relative distance between two terminal subgroups becomes less than a critical radius, which we refer as the capture radius ρ c , by reference to the Smoluchowsky description of first passage reaction kinetics. Fig. 3 shows a bilayer snapshot with its associated terminal subgroups.
Let us now consider a simulation box containing a total of N t lipids, among which N p will be considered a posteriori to be fluorescently labeled, in each leaflet. We denote L the side of a squared simulation box in the x, y directions, with periodic boundary conditions (pbc). We introduce the survival probability P s (ρ c , N p , N t , L; t) as the probability that no excited fluorescent group have come into contact with any of the N p −1 remaining groups in the same leaflet during a time interval t. P s is the outcome of an averaging procedure over both initial conditions (the choice of the fluorescent groups and their actual spatial distribution) and subsequent trajectories (the Brownian displacements of these fluorescent groups). The negative time derivative −dP s /dt represents the probability of "capture" per unit of time, or first passage time distribution. By capture is meant that one of fluorescent group in its ground state penetrates into a sphere of radius ρ c centered around the excited group for the first time during the interval [t, t + dt].
We define in a very similar way P o (ρ c , N p , N t , L; t) the probability that any pair of fluorescent groups taken from randomly selected molecules in two opposite leaflets, have not come closer than a distance ρ c during a time t. The collision dynamics controlling the formation of excimers is order of magnitude faster than the phospholipid flipflop reversal times so that it makes sense to distinguish between the two subpopulations of lipids, and P o effectively controls the influence of interleaflet monomer quenching.
The survival probabilities P s , P o can be efficiently sampled from regular MD trajectories. These dynamical quantities are well defined and related to the molecular displacements properties of the corresponding groups of beads. As we deal with realistic trajectories, we can address situations with transient dynamics, or non Gaussian self-intermediate scattering functions, and possibly improve upon most basic models of molecular diffusion. It remains, however, to relate these probabilities to the excimer/monomer fluorescence ratios which are experimentally obtained.

D. Size-scaling and independent pairs assumptions
Let us consider a lipid bilayer containing a total of N t lipids in its upper leaflet, including a number N p of fluorescent probes. For simplicity, we assume that both leaflets share the same composition, making the bilayer symmetric. The size of the bilayer is supposed to be unaltered by the presence of the probes, so that the size L of the system is directly related to the number of lipids N t by means of the area per lipid a 0 ( N t = L 2 a 0 ).
Under normal illumination conditions, the fraction of time spent in the excited state is tiny, and at any given time the chances of finding two excited monomers in the same portion of membrane is statistically very low. One can therefore consider that one fluorophore at most is excited at a time, and that N p −1 monomers are available to combine as an excimer.
If the motions of the pyrene monomers are decorrelated, the survival probability P s (ρ c , N p , N t ; t) can be reduced to P s (ρ c , 2, N t ; t) Np−1 . This is because N p − 1 different pairs compete independently, and that the effective collision time is the minimum of N p − 1 independent first passage events. Note that the interaction between pyrene monomers is only possible when at least two of them are situated at the immediate vicinity of an excited probe. The outcome of such a situation involving a compact triplet of probes is difficult to ascertain, either enhancing, or disfavoring the formation of a stable excimer. We therefore refer to independent pairs assumption the possibility of factorizing the survival probability when the concentration x of probes is low enough.
When interleaflet interaction is allowed, the survival probability in the independent pairs assumption framework reads: We now address the issue of the sample size dependence. It sounds quite natural to assume that the survival probability depends on size (L or N t ) only throughout the intensive ratio x = N p /N t . If this holds, the survival probability for collisions between opposite leaflets, for instance, should obey the following scaling relation: In the case of collisions between groups in the same leaflet, we observe that the excited pyrene group, i.e. the target, must be removed from the list of moving groups participating to the excimer quenching statistics. This has consequences when N p is small, and the suggested size scaling form for the survival probability within a given leaflet must be slightly modified to account for the "missing" pyrene group.
We refer to these assumptions as the size scaling assumptions. There are reasons to believe, though, that the size scaling assumptions do not hold uniformly. As a matter of fact, collision time distributions do depend explicitly on sample size in 2d Brownian reaction-diffusion models. For instance, the probability of finding a target in a domain of linear size L involves ln(L) corrections with respect to the naive mean-field result. We claim however that these corrections are unlikely to change very much the short time dependence of the survival probability. As far as pyrene excimer formation is concerned, the contribution of long trajectories wandering on large distances ∼ L away from the capture radius are unlikely to contribute to the short time behavior of the survival probability, which turns out to dominate the excimer formation probability. However, they would give rise to significant contributions at longer times.
Combining independent pairs and size scaling assumptions, we can reduce the collision survival probability of an excited pyrenyl group to the survival probability of a single pair sitting in a bilayer of arbitrary size λL: The resulting survival function P c (t) is the global collision survival probability. What remains to be done is to use a scaling factor λ corresponding to a tractable simulation scheme. In the present work, we use λ 2 N t = 256 to sample P s and P o . In addition, we notice that in experimentally relevant situations, λ 2 N t is significantly smaller than the actual number of lipids in the physical system of interest (e.g. a liposome), and the λ 2 term in the exponent of P s can be safely neglected in eq. (4).

E. From survival probabilities to excimer fluorescence intensity
The excimer formation probability is readily obtained from P c (t). A monomer excited at t = 0 forms an excimer during the time interval [t, t + dt] with probability which is the probability of colliding with a ground state monomer while still in the excited state, with τ M standing for the monomer fluorescence lifetime in the ultradilute case (x → 0). The excimer formation probability J E is therefore 1 − J E is the Laplace transform of the survival probability. Naturally, the probability of returning to ground state from the excited monomer state is J M = 1 − J E . The probe ratio x dependence of the monomer/excimer fluorescence intensity curves J M (x), J E (x) comes from the x dependence of the survival probability P c (t), such as expressed in eq. (4). Alternatively, the integrand corresponds to the time-resolved emission intensity of the monomer probe.

F. Intermolecular formation rates
In principle, the survival probability P c must be extracted from the microscopic model for any arbitrary value of x, and the normalized emission intensities J E and J M become dependent on the molecular ratio x. The variation of the fluorescence intensities therefore reports on the concentration of monomeric (or dimeric) probes, and behaves as in a standard titration experiment, except that for practical reasons x cannot be changed during the course of the experiment.
One peculiar value x = x ‡ makes the excimer formation probability equal to 1/2: This value turns out to plays a key role in the interpretation of the experimental titration curves J E (x), J M (x) by Vauhkonen et al. [3] .
Let us assume first that the survival probability P c (t) is exponentially decaying as exp(−t/τ c (x)). In this case, x ‡ is nothing but the concentration for which the collision time τ c (x ‡ ) equals the monomer spontaneous decay time τ M , It is common to treat the excimer formation dynamics as an ordinary bimolecular kinetic process M + M → (M M ) with formation rate K(x), and to take the reverse dissociation as negligible. As the monomer is in excess, the situation corresponds to an ideal pseudo-first order reaction kinetics, consistent with an exponential decay of the excited monomer population. The decrease rate of the excited probes reads d For a pure mean-field, constant rate bimolecular mechanism, K(x) is expected to depend linearly on x, enabling us to rewrite the previous expression as The corresponding normalized emission intensity curves then reduce to simple rational functions.
It can be shown that pseudo-first order excimer formation kinetics and exponential survival behavior are two equivalent assumptions. Conversely, any deviation from linear behavior of the inverse emission intensity 1/J M (x) points to a non exponential survival behavior of the excited monomeric probes, or a non constant rate mechanism. Example of diffusion limited models leading to such time dependent rates are discussed, for instance, in [14]. Our analysis does not make any assumption regarding the kinetic rate coefficients, and the survival probability is sampled from MD trajectories.

G. Multiple passage times
The lack of detailed evidence for the molecular excimer formation mechanism makes it necessary to adopt a phenomenological and probabilistic view on the outcome of the collisional dynamics. Even though the excimer formation is expected to be diffusion dominated, there is nevertheless a chance that collisions does not lead to excimer state with perfect yield.
It is possible to introduce a distance dependent excimer formation rate, such as e.g. the distance dependent energy transfer rate of a Förster resonant pair. In this case, the excimer formation probability depends on the integrated time spent by the pair at close distance.
Denoting by r(t) the relative distance of a single pair, the survival probability of an excited monomer associated with a trajectory {r(t)} reads with q denoting a tunable excimer formation rate, ρ c a characteristic capture radius, and χ a positive function decaying from 1 to 0. The generalization of the survival probability is obtained by averaging over the initial positions and subsequent trajectories r(t): with the same notation as in eq. (3). One notices that the first passage survival probability can be obtained as a limit case of a step function χ[v] = 1 for v ≤ 1 and 0 otherwise, combined with q → ∞. The survival probability for an arbitrary x concentration follows from the size-scaling relation. Functions P s , P o defined as in eq. (14) can be conveniently sampled from molecular dynamics trajectories. This approach to multiple passage dynamics is simple to implement, and contains the first passage capture process as a limit case.
Low values of q correspond to low excimer formation rate. A general feature, as we found, is that the result- ing P s (t) becomes more and more exponential as q decreases, with an average survival time increasing as expected. The resulting inverse probability 1/J E is therefore expected to become more and more linear with x as q decreases, giving us a chance of estimating q. This turned out not to be the case in practice.

H. Estimating the capture radius
There is arbitrariness when it comes to giving a numerical value to the reaction, or capture radius ρ c . In the Martini model, the van der Waals radii of the coarsegrained beads are set to 0.47 nm. Closer approach between beads can only be associated to enthalpic repulsive interactions, which are not supposed to intervene in a diffusion-limited association process.
On the other hand, in order to represent as well as possible the motion of pyrene groups, we bunched together the last beads of each hydrophobic chain groups. The generalized coordinates representing the "pyrene" center of masses, are themselves not materialized as beads and their mutual interaction potential is softer, to some extent, than the one of the original beads. Fig. 4 represents the 2d pair distribution function of group belonging to a given leaflet.
It can be seen that steric hindrance between pyrene center of masses becomes significant for values of r smaller than 0.5 nm. We therefore take ρ c = 0.5 nm as our capture radius value.

I. Real systems vs coarse-grained dynamics
A first question arises, as to which extent coarsegrained dynamics is faithful to the trajectories of the real phospholipid molecules that it is intended to model. It is established that most large scale transport properties such as self-diffusion, transverse viscosity, etc are in fine controlled by cis-trans isomerization dynamics of the lipid alkyl chains [18]. These degrees of freedom are obviously missing at the level of the Martini coarse-grained description, and substituted for by an effective "Kremer-Grest" dynamics of beads and springs polymer chains subject to van-der-Waals interaction. This accounts for the largest fraction of the difference seen between the coarse-grained and atomistic molecular kinetics. Martini simulation are sped up by an acceleration factor f commonly taken of the order of 4. As a by-product, CG molecular dynamics cannot, alone, provides quantitative estimates of the lipid transport coefficients.
The cis-trans chain isomerization moves are temperature dependent, activated processes. As a result, the membrane fluidity strongly depends on temperature, showing Arrhenius dependence in the absence of phase transition. The Martini activation energy is lower than the experimental one. The factor f is therefore likely to depend significantly on temperature. On the other hand, it is not necessary to match precisely the temperature of the CG simulation with the experimental system that one aims at reproducing, provided one restrict ourselves to the same structural phase. We may therefore rely on a kind of time-temperature superposition approximation to adjust to experimental data. Note that irrespective to changes occurring with lipid dynamics, the monomer fluorescence life-time is also a decreasing function of temperature, and must be provided as an external experimental input.
We assume in our approach that the acceleration factor f applies in a uniform way, from short chain reorganization time scales up to long range hydrodynamic limit. In other words, if one can determine f based on the intermediate time regime associated with fluorescence life-time and collision induced excimer formation, one should be able to give a quantitative prediction of the physical lipid diffusion constant D as D = D MD /f , where D M D represents the molecular dynamics diffusion coefficient obtained from the mean-squared displacement. This can be assimilated to a matching procedure.
Another difference between real molecules and CG models lies into the repartition of masses within lipids. The Martini model ascribes an identical mass (72 amu) to all beads, while in reality the phosphate headgroup concentrates more mass than the alkyl chains for comparable steric volumes. As dynamical properties depend on masses (unlike thermodynamics) this could cause differences between real and simulated trajectories. In the present work, we assume that such differences can be disregarded.

IV. NUMERICAL RESULTS AND COMPARISON WITH EXPERIMENTS
A. Diffusion properties of the terminal subgroups, and survival probabilities Survival probabilities were sampled by simulating patches of 512 lipids, respectively POPC at 293K and DOPC at 283K, using 4 trajectories of 150 ns (Martini time). Frames were recorded every 6 ps for the collision statistics, and averages over all the possible pairs of molecules were taken in order to sample P s (t, 2, 256) and P o (t, 1, 256). Trajectories were generated with Gromacs-4.6 [19] using a NVT Nose-Hoover thermostat scheme, and the Martini force-field v2 for the lipids.  (1, 256, t), and then the product P c = P o P s . The latter is expected to rule unrestricted excimer formation. Introducing collisions between opposite leaflets increases significantly the survival decay rate. Fig. 6 shows the collisional survival probability P c (2, 256, t) of the chain subgroups, that is central to our reaction kinetics modeling, on a semi-logarithmic plot. As discussed in the previous section, such non exponential behavior can be related to deviation from the Stern-Volmer linear plot of the titration curves. A direct test of the size-scaling approximation is provided in Fig. 7, when survival probabilities obtained from a single pair of monomers embedded respectively in leaflets of 256 and 1024 lipids are displayed. The scaled survival curve P s (2, 1024, t) 4 matches P s (2, 256, t) at short times, while departing from it at longer times. This expected behavior originates from trajectories showing a large separation between reacting pairs, allowed in the large system but forbidden (due to periodic boundary conditions) in the small system. It is also a signature of the size dependence of Brownian diffusion limited reaction dynamics in two-dimensional systems. Clearly, departure from the size scaling assumptions is visible on this curve, at large times. However, the most relevant region, as far as excimer formation is concerned, is the short time regime, and this is especially true at large DOPC POPC  probe concentrations x. Therefore, we consider that the size-scaling assumption is valid in our case. Fig. 8 shows that survival probabilities obtained from trajectories at different temperatures almost superimpose, once a dimensionless time variable u = Dt/a 0 is used on the horizontal axis.

B. Experimental titration curves
The experimental data on monomer to excimer emission ratio are taken from Fig. 3 in reference [3] (read from the graph) and displayed in Table I. The experimental values of the monomer fluorescence lifetime, which depends on temperature, are taken from Table 1 in reference [4]. The lifetimes of the isolated monomers are estimated from Fig. 5 in [3] to be respectively 140 ns for POPC at 20 • C and 160 ns for DOPC at 10 • C. It must be stressed that the final results are very much dependent on the actual τ M values.

Lipid
Original D MD D MD D' POPC (293 K) 11 1.8 4.0 DOPC (283 K) 6 1.0 1. 6   TABLE II. Comparison of diffusion constants, obtained with the milling crowd analysis of [3] (column 2) and with the current approach, assuming reaction at first contact and allowing for interleaflet association (column 3) or excluding interleaflet association (column 4).
C. Fit of the excimer formation dynamics Fig. 9(top left) shows the agreement between the experimental values and the numerical prediction extracted from the survival probability of Fig. 6. The best agreement is obtained with "numerical" fluorescence lifetime τ M,CG equals to 10.7 ns. The ratio between the experimental and the numerical lifetimes determines the acceleration factor f = τ M /τ M,CG = 15, given that τ M = 160 ns. The coarse-grained diffusion constant D CG is assumed to scale as f D, with f the acceleration factor and D the real time diffusion constant. This approach predicts a value D = 1.6 µm.s −1 for DOPC molecules in a bilayer at 10 • C, if no interleaflet association are allowed.
When adjustment is made with POPC trajectories at 20 K, as shown in Fig. 9(top right), the best fit is obtained with a numerical lifetime τ M,CG equal to 15.6 ns. The acceleration factor is now f = 9 and the associated diffusion constant found for POPC at 20 • C is D = 4.0 µm.s −1 . Fig. 9(bottom left) shows the best agreement between the experimental and numerical curves, when the survival probability P c (t) = P s (t) × P o (t) corresponding to unrestricted excimer formation is used. As the decay is now faster than in the previous case, a larger acceleration factor must be used to match the MD trajectories to the observed excimer formation rate. We obtain for DOPC at 10 • C a factor α = 26, leading to a diffusion constant D = 1.0 µm.s −1 . A similar adjustment to the experimental data leads, for a POPC bilayer at 20 • C, to an acceleration factor f = 20 and D = 1.8 µm.s −1 , shown in Fig. 9(bottom right).

V. DISCUSSION
We first observe that our procedure reproduce to some extend the upward curvature seen on the 1/J M (x) plots, even though a larger curvature would be required to optimally fit the data.
In Table II should be compared with our column 4. We find a significant reduction of the diffusion constant compared with their analysis.
Our diffusion constant in the case of allowed interleaflet formation is even smaller. On the one hand, if one trusts the coarse-grained MD trajectories, there is no reason to discard interleaflet association which occurs at a non negligible rate. In other words, P c (t) differs enough from P s (t) to alter significantly the final results. The pyrene probes by Vauhkonen et al. are the same, irrespective of the bilayers under consideration. Interleaflet association is expected to occur significantly for short lipids (e.g. DMPC) but to be disfavored for longer lipid (e.g. DOPC, DPPC, POPC) if there is a barrier preventing pyrene groups from interdigitating. Our coarse-grained model does not show evidence to this, but there is the possibility that it misses this point. United atoms simulations of pyrenyl groups grafted to saturated lipid chains in a DPPC fluid environment, as well as dipyrenyl groups grafted at the end of chains with variable length have been published [20][21][22]. These studies suggests that interdigitation of the pyrenyl group is possible if grafted to a 10 carbons length chain, as in the experiments considered here, and unlikely for shorter chains. Note that results obtained for saturated lipid chain may not be transposable to POPC and DOPC membranes.
On the other hand, the results obtained by restricting pyrene association to same leaflet are closer from the body of published experimental coefficients. The Handbook of phospholipid bilayers published by Marsh [23] reports for DOPC at 25 • C the following two values: 6.3µm 2 .s −1 (fluorescence correlation spectroscopy) and 1.8 µm 2 .s −1 (electron spin resonance). Using an energy activation E a ∼ 11 kJ.mol −1 , these values transpose to respectively 5 and 1.4 µm 2 .s −1 . It appears that the experimental diffusion coefficients reported in [23] are widely scattered, and much higher values of D are even reported. For POPC at 20 • C a FRAP value of 3.4 µm 2 .s −1 is given. Again, it is consistent with our value in the absence of interleaflet association. Moreover, our new value is closer to FRAP than the original interpretation of the same data.
To conclude with diffusion coefficients, we find that when interleaflet association is allowed, our coefficients are a factor 2 smaller than those corresponding to fluorescent probe diffusion (FRAP or FCS). A better agreement is obtained by restricting pyrene association to probes occupying the same leaflet.
An outcome of the present work is the determination of the acceleration factor f of a Martini CG lipid bilayer model. At 10 • C, f lies between 15 and 25, depending on interleaflet association. At 20 • C it is reduced to 10-20. This is significantly higher than f = 4 which is sometimes assumed. The f = 4 factor corresponds indeed to the equivalent solvent (water) diffusion model, but seems higher as far as lipid motion is concerned. As the CG activation energy (14 kJ/mol) is likely to differ from the real activation energy, the agreement between CG model and reality should improve in the temperature range 300-330 K for which this model is supposed to perform best. Working at low temperatures enlarge the dynamical spread between numerical and real trajectories.
Coarse-graining raises the question of whether one can trust the resulting time correlation functions. There is no reason to assume that the acceleration factor applies uniformly for each time interval. If one denotes C CG AB (t) = A(t)B(0) the correlation function of a pair of observables, and C AB (t) = A(t)B(0) the experimental (or all-atom simulated) counterpart, a correspondence C CG AB (T (t)) = C AB (t) is expected, with possibly a nonlinear monotonous time correspondence T (t) between the two correlations.
An homogeneous acceleration factor corresponds to T (t) = t/f . An inhomogeneous correspondence could explain why is the short time collisional dynamics faster than the long time diffusion process. Our work sets bounds on such a phenomenon. If interleaflet association is forbidden, the agreement between numerical and experimental values of D is consistent with an homogeneous acceleration factor. If interleaflet association is allowed, the time-scale dependent effective factor f = 25 (DOPC) or f = 20 (POPC) that acts on short time separation, reduced by a factor 2, would be again consistent with a diffusion coefficient of respectively 2 and 3.6 µm 2 .s −1 .
Another way of slowing down the collision dynamics would be to challenge the diffusion limited character of the excimer formation, by requiring repeated collisions between monomers prior to complexation. This effect can be introduced by changing the excimer formation rate q in equation (14). For DOPC, interleaflet association allowed, a diffusion coefficient D = 2.2 µm 2 .s −1 can be obtained by reducing q to a value such that it would take 1.1 ns of real time for the excimer to form, with both monomers maintained at close vicinity. This extra time could be related to a need of the terminal subgroups to align parallel, prior to excimer formation.
A drawback of slowing back in such a way the dynamics is that the effective survival probability becomes closer and closer to exponential shape, and therefore does not lead to the desired upward curvature of the 1/J M (x) plot. Therefore, we do not find that it improves the agreement with experimental data to introduce such an intrinsic excimer formation time. A clue to whether such a delay is necessary, and its order of magnitude could be obtained from time-resolved excimer formation.

VI. CONCLUSION
We have presented a very general and efficient method to derive association kinetic coefficients in the case of diffusion limited processes, based on realistic molecular motions in lipid bilayers. The time-dependent reaction rate is derived from a survival probability, which itself comes from a posteriori trajectory analysis. A virtual relabeling of the simulated molecules enables us to efficiently gather collision statistics. No assumptions are made regarding the Brownian nature of the molecular motion, and intramolecular fluctuation effects are taken into account, within the limits of the coarse-graining procedure.
Our main assumptions are that lipid probes behave similarly to the major lipid components, and do not interact strongly with each other. We also require the dynamical acceleration factor of the coarse-grained numerical model to be uniform over the relevant time scales.
It is first shown that simulations at different sizes and temperatures can be compared, and used for the purpose of studying the collision dynamics between virtually labeled lipids. A likely capture radius ρ c = 0.5 nm is obtained from a pair correlation analysis of the reactive groups.
Two main situations are considered, depending on whether interleaflet excimer association is allowed or not. Based on CG simulations considerations, interleaflet association should be allowed, but lead to small values of the diffusion coefficients. If one forbids interleaflet association, given that there could barriers opposing it that our CG model misses, the resulting diffusion constants fall in a range comparable to other independent experimental procedures. Between these extremes, there is also the possibility that interleaflet excimer formation occurs, but comparatively less frequently than intraleaflet excimers.
In any case, our approach leads to a significant reduc-tion of the diffusion coefficients compared to those originally published. We suggest that the coefficients found in similar diffusion limited processes, which often lie above the values obtained from long range fluorescent probes diffusion, should be reanalyzed along the lines of our present approach. We suspect that these diffusion coefficients are overestimated due to analysis bias, as our example shows. Our standing is that we have possibly explained the discrepancy between these two different approaches to lipid diffusion determination. More modeling work is needed in order to determine the degree of interdigitation of pyrene groups pertaining to different leaflets, as well as regarding the effect of non trivial dispersion forces between resonant fluorescent groups.
The application of this technique to complex membranes with dynamical heterogeneities would certainly be interesting, now that a convenient framework for diffusion limited reaction rate analysis, based on molecular dynamics, is available.
In a mean-field, constant rate bimolecular mechanism, K(x) is proportional to x and therefore again It is clear from eqs. (A3) and (A6) that both approaches coincide.

Multiple passage times
The possible strategy, which we did not consider in our main paper, consists in introducing a generalization P (n) C of the survival probability P C . Let P (n) C be the probability that the excited monomer has suffered exactly n collisions during the excitation time t. The first term of the sequence P If trajectories and excimer formation reaction are uncorrelated, one can introduce an effective survival probability with W E the probability to form an excimer at a given encounter, such that eq. (13) applies with P eff (t) substituted for P C . When using this strategy in the Poissonian case, one finds that the effective survival probability remains exponential with an effective decay time (1 − W E )τ c . As a byproduct, the titration curve is strictly independent of W E . Therefore it is not possible to determine W E in the case of an uncorrelated, Poissonian collision dynamics process. Irrespective of W E the inverse monomer probability reads 1/J M (x) = 1 + x/x ‡ . Some deviation from pure exponential behavior of the survival function is necessary if one wants to observe and fit a nonlinear curve 1/J M (x).