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

P. Ayoub, F. Thalmann∗ Institut Charles Sadron, Université de Strasbourg, CNRS
UPR 22, 23 rue du Loess, Strasbourg Cedex, F-67037, France.
<EMAIL_ADDRESS>
###### Abstract

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
1974_Galla_Sackmann , 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 2007_Marrink_deVries . 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. 1990_Vauhkonen_Eisinger ;
1990_Sassaroli_Eisinger . 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 1975_Birks . 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 membranes. 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 2002_Somerharju . A larger
class of diffusion-limited excitation or deexcitation processes in membranes
is covered in a review by Melo and Martins 2006_Melo_Martins .

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
$\tau_{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
1974_Galla_Sackmann ; 1990_Vauhkonen_Eisinger ; 1990_Sassaroli_Eisinger ;
1979_Galla_Sackmann ; 1996_Martins_Melo ; 2001_Martins_Melo ;
2001_Novikov_Visser . 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 1974_Naqvi ;
1974_Naqvi_Chapman ; 1983_Torney_McConnell ; 1989_Szabo . 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 coarse-grained model for lipid
bilayers, the Martini force field, and to analyze the excimer fluorimetric
data based on realistic molecular dynamics (MD) trajectories.

##  III Theory

###  III.1 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
2010_Zhou .

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 Smoluchovsky,
considers two spheres following ideal Brownian trajectories and reacting as
soon as their relative distance drops below a critical value: the capture
radius $\rho_{c}$. Naqvi derived a Smoluchowsky rate in the two dimensional
case 1974_Naqvi . Both approaches predict a time-dependent 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
1983_Torney_McConnell .

The experimental results reported in 1990_Sassaroli_Eisinger  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 mean-squared 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
coarse-graining 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.

###  III.2 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 relatively 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 1975_Birks ; 1964_Stephen . 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 $\to$ diffusion $\to$
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.

###  III.3 Survival probabilities

Figure 1: Example of coarse-grained lipid molecule (DOPC) with bead labels,
and two terminal subgroups during a “collision event”.

Figure 2: Detail of two DOPC molecules, each in one leaflet, with terminal
subgroups superimposed on top of them.

Figure 3: Side view of a DOPC bilayer, with superimposed terminal subgroups
from the upper and lower leaflets (using two different colors depending on the
leaflet to which the molecule belongs).

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 $\rho_{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}(\rho_{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 $-\mathrm{d}P_{s}/\mathrm{d}t$
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 $\rho_{c}$ centered around the
excited group for the first time during the interval $[t,t+\mathrm{d}t]$.

We define in a very similar way $P_{o}(\rho_{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
$\rho_{c}$ during a time $t$. The collision dynamics controlling the formation
of excimers is order of magnitude faster than the phospholipid flip-flop
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.

###  III.4 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}(\rho_{c},N_{p},N_{t};t)$ can be reduced to
$P_{s}(\rho_{c},2,N_{t};t)^{N_{p}-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:

$P_{\mathrm{c}}(t)=P_{s}(\rho_{c},2,N_{t},L;t)^{N_{p}-1}P_{o}(\rho_{c},1,N_{t},L;t)^{N_{p}}.$
(1)

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:

$\displaystyle P_{o}(\rho_{c},N_{p},N_{t};t)$ $\displaystyle=$ $\displaystyle
P_{o}(\rho_{c},\lambda^{2}N_{p},\lambda^{2}N_{t};t)$ (2) $\displaystyle=$
$\displaystyle P_{o}(\rho_{c},\lambda^{2}xN_{t},\lambda^{2}N_{t};t).$

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.

$\displaystyle P_{s}(\rho_{c},N_{p},N_{t};t)$ $\displaystyle=$ $\displaystyle
P_{s}(\rho_{c},1+\lambda^{2}(N_{p}-1),\lambda^{2}N_{t};t)$ $\displaystyle=$
$\displaystyle P_{s}(\rho_{c},1+\lambda^{2}(xN_{t}-1),\lambda^{2}N_{t};t).$

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 $\sim 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 $\lambda
L$:

$\displaystyle P_{\mathrm{c}}(t)$ $\displaystyle=$ $\displaystyle
P_{s}(\rho_{c},2,\lambda^{2}N_{t};t)^{\lambda^{2}xN_{t}-\lambda^{2}}$ (4)
$\displaystyle\;\;\times\;P_{o}(\rho_{c},1,\lambda^{2}N_{t};t)^{\lambda^{2}xN_{t}},$
$\displaystyle\simeq$ $\displaystyle
P_{s}(\rho_{c},2,\lambda^{2}N_{t};t)^{\lambda^{2}xN_{t}}$
$\displaystyle\;\;\times\;P_{o}(\rho_{c},1,\lambda^{2}N_{t};t)^{\lambda^{2}xN_{t}}.$

The resulting survival function $P_{\mathrm{c}}(t)$ is the global collision
survival probability.

What remains to be done is to use a scaling factor $\lambda$ corresponding to
a tractable simulation scheme. In the present work, we use
$\lambda^{2}N_{t}=256$ to sample $P_{s}$ and $P_{o}$. In addition, we notice
that in experimentally relevant situations, $\lambda^{2}N_{t}$ is
significantly smaller than the actual number of lipids in the physical system
of interest (e.g. a liposome), and the $\lambda^{2}$ term in the exponent of
$P_{s}$ can be safely neglected in eq. (4).

###  III.5 From survival probabilities to excimer fluorescence intensity

The excimer formation probability is readily obtained from
$P_{\mathrm{c}}(t)$. A monomer excited at $t=0$ forms an excimer during the
time interval $[t,t+\mathrm{d}t]$ with probability

$-\frac{\mathrm{d}P_{\mathrm{c}}}{\mathrm{d}t}\exp\left(-\frac{t}{\tau_{M}}\right),$
(5)

which is the probability of colliding with a ground state monomer while still
in the excited state, with $\tau_{M}$ standing for the monomer fluorescence
lifetime in the ultra-dilute case ($x\to 0$). The excimer formation
probability $J_{E}$ is therefore

$\displaystyle J_{E}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}-\frac{\mathrm{d}P_{\mathrm{c}}}{\mathrm{d}t}\exp\left(-\frac{t}{\tau_{M}}\right)\mathrm{d}t$
(6) $\displaystyle=$ $\displaystyle
1-\int_{0}^{\infty}P_{\mathrm{c}}(t)\exp\left(-\frac{t}{\tau_{M}}\right)\frac{\mathrm{d}t}{\tau_{M}}.$

$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_{\mathrm{c}}(t)$, such as expressed
in eq. (4). Alternatively, the integrand

$\frac{1}{\tau_{M}}P_{\mathrm{c}}(t)\exp\left(-\frac{t}{\tau_{M}}\right)$ (7)

corresponds to the time-resolved emission intensity of the monomer probe.

###  III.6 Intermolecular formation rates

In principle, the survival probability $P_{\mathrm{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^{\ddagger}$ makes the excimer formation probability
equal to $1/2$:

$J_{M}(x^{\ddagger})=J_{E}(x^{\ddagger})=1/2.$ (8)

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.
1990_Vauhkonen_Eisinger  .

Let us assume first that the survival probability $P_{\mathrm{c}}(t)$ is
exponentially decaying as $\exp(-t/\tau_{c}(x))$. In this case, $x^{\ddagger}$
is nothing but the concentration for which the collision time
$\tau_{c}(x^{\ddagger})$ equals the monomer spontaneous decay time $\tau_{M}$,

$\tau_{c}(x^{\ddagger})=\tau_{M}.$ (9)

It is common to treat the excimer formation dynamics as an ordinary
bimolecular kinetic process $M+M^{\star}\to(MM)^{\star}$ 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
$\mathrm{d}[M^{\star}]=-(K(x)+\tau_{M}^{-1})[M^{\star}]\mathrm{d}t$, leading
to the monomer emission intensity

$\frac{1}{J_{M}(x)}=1+\tau_{M}K(x).$ (10)

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

$\frac{1}{J_{M}(x)}=1+\frac{x}{x^{\ddagger}}.$ (11)

The corresponding normalized emission intensity curves then reduce to simple
rational functions.

$J_{M}(x)=\frac{x^{\ddagger}}{x+x^{\ddagger}}\;;\;J_{E}(x)=\frac{x}{x+x^{\ddagger}}.$
(12)

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
1983_Torney_McConnell . Our analysis does not make any assumption regarding
the kinetic rate coefficients, and the survival probability is sampled from MD
trajectories.

###  III.7 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

$\exp\left(-q\int_{0}^{t}\chi\left[\frac{r(t^{\prime})}{\rho_{c}}\right]\,\mathrm{d}t^{\prime}\right),$
(13)

with $q$ denoting a tunable excimer formation rate, $\rho_{c}$ a
characteristic capture radius, and $\chi$ 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)$:

$P_{s}(q,\rho_{c},2,\lambda^{2}N_{t};t)=\\\
\left\langle\exp\left(-q\int_{0}^{t}\chi\left[\frac{r(t^{\prime})}{\rho_{c}}\right]\,\mathrm{d}t^{\prime}\right)\right\rangle_{\\{r(t)\\}},$
(14)

with the same notation as in eq. (LABEL:eq:survivalProba). One notices that
the first passage survival probability can be obtained as a limit case of a
step function $\chi[v]=1$ for $v\leq 1$ and $0$ otherwise, combined with
$q\to\infty$. 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 resulting $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.

###  III.8 Estimating the capture radius

There is arbitrariness when it comes to giving a numerical value to the
reaction, or capture radius $\rho_{c}$. In the Martini model, the van der
Waals radii of the coarse-grained 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.

Figure 4: Pair distribution of terminal subgroups. The vertical dashed line
indicates the horizontal separation $r=0.5$.

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
$\rho_{c}=0.5$ nm as our capture radius value.

###  III.9 Real systems vs coarse-grained dynamics

A first question arises, as to which extent coarse-grained 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
Cevc_Marsh_PhospholipidBilayers . 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_{\mathrm{MD}}/f$, where $D_{MD}$ 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

###  IV.1 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 2008_Hess_Lindahl  using a NVT Nose-Hoover thermostat scheme, and
the Martini force-field v2 for the lipids.

Figure 5: Survival probabilities restricted to same leaflet, opposite leaflet
and product of both. Trajectories are from DOPC at 283 K, with a capture
radius $\rho_{c}=0.5$ nm and frames taken every 6 ps of simulation time.

Fig. 5 shows the survival probability when collisions are restricted to
molecules belonging to the same leaflet $P_{s}(2,256,t)$, to opposite leaflets
$P_{o}(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.

Figure 6: Survival probability $P_{\mathrm{c}}(t)$ when collisions are
unrestricted (DOPC at 283 K). The vertical axis is logarithmic, and a straight
line representing the asymptotic behavior $0.926\exp(-0.00196t)$ ($t$ in ns)
is also shown. The straight line goes from 0.926 ($t$=0) to 0.687 ($t$=150
ns), and the distance to upper curve represents a deviation from pure
exponential behavior of $P_{\mathrm{c}}$. The curve represented above
corresponds to a concentration $x=1/256$, and the non-exponential behavior at
short times is reinforced as the concentration $x$ increases.

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.

Figure 7: Size-scaling approximation for DOPC molecules, using 512 and 2048
lipids in the simulations.

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
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.

Figure 8: Time temperature superposition of the survival probability for DOPC
bilayers. Simulated temperatures are 283, 300 and 330K.

###  IV.2 Experimental titration curves

The experimental data on monomer to excimer emission ratio are taken from Fig.
3 in reference 1990_Vauhkonen_Eisinger  (read from the graph) and displayed in
Table 1. The experimental values of the monomer fluorescence lifetime, which
depends on temperature, are taken from Table 1 in reference
1990_Sassaroli_Eisinger . The lifetimes of the isolated monomers are estimated
from Fig. 5 in 1990_Vauhkonen_Eisinger  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 $\tau_{M}$ values.

| DOPC |  |  | POPC |
---|---|---|---|---|---
$x$ | $1/J_{M}(x)$ | $J_{M}(x)$ | $x$ | $1/J_{M}(x)$ | $J_{M}(x)$
0.003 | 1.02 | 0.98 | 0.003 | 1.03 | 0.97
0.012 | 1.09 | 0.92 | 0.012 | 1.14 | 0.88
0.034 | 1.31 | 0.76 | 0.033 | 1.52 | 0.66
0.055 | 1.58 | 0.63 | 0.056 | 1.93 | 0.52
0.102 | 2.45 | 0.41 | 0.102 | 3.25 | 0.31
Table 1: Monomer fraction $J_{M}(x)$ vs probe concentration $x$

###  IV.3 Fit of the excimer formation dynamics

|
---|---
|
Figure 9:  Top left: fit of $1/J_{M}(x)$ in the DOPC case when association is
restricted to molecules within a single leaflet. Top right: fit of
$1/J_{M}(x)$ in the POPC case when association is restricted to molecules
within a single leaflet. Bottom left: fit of $1/J_{M}(x)$ in the DOPC case
when association originates from both leaflets. Bottom right: fit of
$1/J_{M}(x)$ in the POPC case when association originates from both leaflets.

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
$\tau_{M,CG}$ equals to 10.7 ns. The ratio between the experimental and the
numerical lifetimes determines the acceleration factor
$f=\tau_{M}/\tau_{M,CG}=15$, given that $\tau_{M}=160$ ns. The coarse-grained
diffusion constant $D_{CG}$ is assumed to scale as $fD$, with $f$ the
acceleration factor and $D$ the real time diffusion constant. This approach
predicts a value $D=1.6~{}\mu\mathrm{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 $\tau_{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~{}\mu\mathrm{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)\times
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 $\alpha=26$, leading to a diffusion
constant $D=1.0~{}\mu\mathrm{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~{}\mu\mathrm{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.

Lipid | Original $D$ | MD $D$ | MD $D$’
---|---|---|---
POPC (293 K) | 11 | 1.8 | 4.0
DOPC (283 K) | 6 | 1.0 | 1.6
Table 2: Comparison of diffusion constants, obtained with the milling crowd
analysis of 1990_Vauhkonen_Eisinger  (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).

In Table 2, we compare our new values to those originally published by
Vauhkonen et al. (Table 1 of ref. 1990_Vauhkonen_Eisinger ). Strictly
speaking, the modeling of Vauhkonen et al excludes interleaflet excimer
formation, and their values 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
2006_Repakova_Vattulainen ; 2014_Franova_Ollila ; 2014_Franova_Vattulainen .
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
Marsh_HandbookLipidBilayers2  reports for DOPC at 25∘C the following two
values: 6.3$\mu\mathrm{m}^{2}.s^{-1}$ (fluorescence correlation spectroscopy)
and 1.8 $\mu\mathrm{m}^{2}.s^{-1}$ (electron spin resonance). Using an energy
activation $E_{a}\sim 11~{}\mathrm{kJ.mol}^{-1}$, these values transpose to
respectively 5 and 1.4 $\mu\mathrm{m}^{2}.s^{-1}$. It appears that the
experimental diffusion coefficients reported in Marsh_HandbookLipidBilayers2
are widely scattered, and much higher values of $D$ are even reported. For
POPC at $20^{\circ}$C a FRAP value of 3.4 $\mu\mathrm{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^{\circ}$C, $f$ lies between 15
and 25, depending on interleaflet association. At $20^{\circ}$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_{AB}^{CG}(t)=\langle A(t)B(0)\rangle$ the correlation function of a pair of
observables, and $C_{AB}(t)=\langle A(t)B(0)\rangle$ the experimental (or all-
atom simulated) counterpart, a correspondence $C_{AB}^{CG}(T(t))=C_{AB}(t)$ is
expected, with possibly a non-linear 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 $\mu\mathrm{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~{}\mu\mathrm{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 $\rho_{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 reduction 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.

## Author contributions

P.A. did the simulation and numerical work, the data fitting and composed the
figures (4-9). F.T. wrote the manuscript. Both authors contributed to the
theoretical model and to figures (1-3).

## Acknowledgments

F.T. would like to thank G. Duportail, C.M. Marques and J. Léonard for
discussions. The authors acknowledge support from the High Performance Cluster
Equip@Meso of the University of Strasbourg for providing computational
resources.

## References

  * (1) Hans-Joachim Galla and Erich Sackmann.  Lateral diffusion in the hydrophobic region of membranes: use of pyrene excimers as optical probes.  Biochimica et Biophysica Acta (BBA) - Biomembranes, 339(1):103 – 115, 1974.
  * (2) Siewert J. Marrink, H. Jelger Risselada, Serge Yefimov, D. Peter Tieleman, and Alex H. de Vries.  The martini force field: a coarse grained model for biomolecular simulations.  The Journal of Physical Chemistry B, 111(27):7812–7824, 2007.  PMID: 17569554.
  * (3) M. Vauhkonen, M. Sassaroli, P. Somerharju, and J. Eisinger.  Dipyrenylphosphatidylcholines as membrane fluidity probes. relationship between intramolecular and intermolecular excimer formation rates.  Biophysical Journal, 57(2):291 – 300, 1990.
  * (4) M. Sassaroli, M. Vauhkonen, D. Perry, and J. Eisinger.  Lateral diffusivity of lipid analogue excimeric probes in dimyristoylphosphatidylcholine bilayers.  Biophysical Journal, 57(2):281 – 290, 1990.
  * (5) J. B. Birks.  Excimers.  Reports on Progress in Physics, 38(8):903, 1975.
  * (6) Pentti Somerharju.  Pyrene-labeled lipids as tools in membrane biophysics and cell biology.  Chemistry and Physics of Lipids, 116(1-2):57 – 74, 2002.
  * (7) Eurico Melo and Jorge Martins.  Kinetics of bimolecular reactions in model bilayers and biological membranes. a critical review.  Biophysical Chemistry, 123(2-3):77 – 94, 2006.
  * (8) H.-J. Galla, W. Hartmann, U. Theilen, and E. Sackmann.  On two-dimensional passive random walk in lipid bilayers and fluid pathways in biomembranes.  The Journal of Membrane Biology, 48:215–236, 1979.
  * (9) Jorge Martins, Winchil L. C. Vaz, and Eurico Melo.  Long-range diffusion coefficients in two-dimensional fluid media measured by the pyrene excimer reaction.  The Journal of Physical Chemistry, 100(5):1889–1895, 1996.
  * (10) Jorge Martins and Eurico Melo.  Molecular mechanism of lateral diffusion of py10-pc and free pyrene in fluid {DMPC} bilayers.  Biophysical Journal, 80(2):832 – 840, 2001.
  * (11) EG Novikov and AJWG Visser.  Inter- and intramolecular dynamics of pyrenyl lipids in bilayer membranes from time-resolved fluorescence spectroscopy.  Journal of Fluorescence, 11(4):297–305, 2001.
  * (12) K. Razi Naqvi.  Diffusion-controlled reactions in two-dimensional fluids: discussion of measurements of lateral diffusion of lipids in biological membranes.  Chemical Physics Letters, 28(2):280 – 284, 1974.
  * (13) K.Razi Naqvi, Jean-Paul Behr, and D. Chapman.  Methods for probing lateral diffusion of membrane components: triplet-triplet annihilation and triplet-triplet energy transfer.  Chemical Physics Letters, 26(3):440 – 444, 1974.
  * (14) D. C. Torney and H. M. McConnell.  Diffusion-limited reaction rate theory for two-dimensional systems.  Proceedings of the Royal Society of London. A, 387(1792):147–170, 1983.
  * (15) Attila Szabo.  Theory of diffusion-influenced fluorescence quenching.  The Journal of Physical Chemistry, 93(19):6929–6939, 1989.
  * (16) Huan-Xiang Zhou.  Rate theories for biologists.  Quarterly Review of Biophysics, 43(2):219–293, 2010.
  * (17) M. J. Stephen.  First-order dispersion forces.  The Journal of Chemical Physics, 40(3):669–673, 1964.
  * (18) Gregor Cevc and Derek Marsh.  Phospholipid Bilayers. Physical Principles and Models.  John Wiley & Sons, New-York, 1987.
  * (19) Berk Hess, Carsten Kutzner, David van der Spoel, and Erik Lindahl.  Gromacs 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation.  Journal of Chemical Theory and Computation, 4(3):435–447, 2008\.
  * (20) J. RepÃ¡kovÃ¡, J. M. Holopainen, M. Karttunen and I. Vattulainen.  Influence of Pyrene-Labeling on Fluid Lipid Membranes  Journal of Physical Cheistry B, 110:15403-15410, 2006.
  * (21) M. D. FraÅ^ovÃ¡, I. Vattulainen and O. S. Ollila.  Can pyrene probes be used to measure lateral pressure profiles of lipid membranes? Perspective through atomistic simulations  BBA-Biomembranes, 1838:1406-1411, 2014.
  * (22) M. D. FraÅ^ovÃ¡, and J. RepÃ¡kovÃ¡, and J. M. Holopainen and I. Vattulainen.  How to link pyrene to its host lipid to minimize the extent of membrane perturbations and to optimize pyrene dimer formation  Chemistry and Physics of Lipids, 177:19-25, 2014.
  * (23) Derek Marsh.  Handbook of Lipid Bilayers.  CRC Press, Boca Raton, 2nd edition, 2013.


# An Off-Axis Relativistic Jet Model for the Type Ic supernova SN 2007gr

M. Xu11affiliation: Department of Astronomy, Nanjing University, Nanjing
210093, China<EMAIL_ADDRESS>22affiliation: Yukawa Institute for Theoretical
Physics, Oiwake-cho, Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan
33affiliation: Key Laboratory of Modern Astronomy and Astrophysics (Nanjing
University), Ministry of Education, China , S. Nagataki22affiliation: Yukawa
Institute for Theoretical Physics, Oiwake-cho, Kitashirakawa, Sakyo-ku, Kyoto
606-8502, Japan , and Y. F. Huang11affiliation: Department of Astronomy,
Nanjing University, Nanjing 210093, China<EMAIL_ADDRESS>33affiliation: Key
Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry
of Education, China

###### Abstract

We propose an off-axis relativistic jet model for the Type Ic supernova SN
2007gr. Most of the energy ($\sim 2\times 10^{51}$ erg) in the explosion is
contained in non-relativistic ejecta which produces the supernova. The optical
emission is coming from the decay process of $\rm{}^{56}Ni$ synthesized in the
bulk SN ejecta. Only very little energy ($\sim 10^{48}$ erg) is contained in
the relativistic jet with initial velocity about $0.94$ times the speed of
light. The radio and X-ray emission comes from this relativistic jet. With
some typical parameters of a Wolf-Rayet star (progenitor of Type Ic SN), i.e.,
the mass loss rate $\dot{M}=1.0\times 10^{-5}~{}M_{\odot}~{}\rm yr^{-1}$ and
the wind velocity $v_{\rm w}=1.5\times 10^{3}~{}\rm km~{}s^{-1}$ together with
an observing angle of $\theta_{\rm obs}=63.3^{\circ}$, we can obtain the
multiband light curves that fit the observations well. All the observed data
are consistent with our model. Thus we conclude that SN 2007gr contains a weak
relativistic jet and we are observing the jet from off-axis.

gamma-ray bursts: general - supernovae: individual (SN 2007gr)

##  1 Introduction

Gamma-ray bursts (GRBs) are intense flashes of gamma-ray radiation in the
Universe (for recent reviews, see: Zhang 2007; Gehrels et al. 2009). It is
widely believed that outflows of GRBs are accelerated to ultra-relativistic
speeds (Mészáros 2002) and usually collimated with small jet anlges (Frail et
al. 2001). On the other side, core-collapse supernovae (SNe) are the explosive
deaths of massive stars that occur when their iron cores collapse to form
neutron stars or black holes (Wilson 1971; Barkat et al. 1974; Wheeler &
Levreault 1985; Woosley & Janka 2005; Woosley & Bloom 2006; Nagataki et al.
2007; Nagataki 2009; Nagataki 2010). According to their spectra, core-collapse
SNe are classified as Type Ib, Ic or Type II SNe (Wheeler 1993; Wheeler et al.
1993; Filippenko 1997). SN explosions will release massive ejecta and usually
they are isotropic and non-relativistic.

Recent observations have revealed that many nearby GRBs are associated with
core-collapse SNe. Examples of such association includes GRB 980425/SN 1998bw
(Galama et al. 1998), GRB 030329/SN 2003dh (Berger et al. 2003), GRB 031203/SN
2003lw (Cobb et al. 2004), GRB 060618/SN 2006aj (Campana et al. 2006), GRB
091127/SN 2009nz (Cobb et al. 2010), and GRB 100316D/SN 2010bh (Fan et al.
2010) etc. These GRBs are usually soft in $\gamma$-ray spectra and are
ubiquitously longer than $2$ seconds. Thus they belong to the so called
long/soft GRBs.

The observed GRB-connected SNe are all Type Ic SNe. The most favored
progenitors for Type Ic supernovae are Wolf-Rayet stars (Maeder & Lequeux
1982; Begelman & Sarazin 1986; Woosley & Bloom 2006). However, the GRB-
connected Type Ic SNe should be different from ordinary Type Ic SNe (Soderberg
et al. 2006), because they need to launch relativistic jets to produce the
bursts of $\gamma$-rays. The kinetic energy of these SNe appears to be greater
than that of ordinary SNe. In some SNe associated with GRBs, most of the
explosion energy is in non-relativistic ejecta which produces the supernova,
while only little energy is in the relativistic jets which are responsible for
making GRBs and their afterglows (Woosley & Bloom 2006).

While GRB-connected SN explosions can produce relativistic jets with Lorentz
factor as large as $\sim 100$, we believe that there were still some Type Ic
SNe that could only produce midly relativistic jets with initial Lorentz
factor of a few (Huang et al. 2002; Granot & Loeb 2003). It has been argued
that this kind of low Lorentz factor jets will produce UV or soft X-ray
transients but not GRBs (Huang et al. 2002; Xu et al. 2010). The interesting
Type Ic supernova SN 2009bb, which is identified as a relativistic ejecta
without a detected GRB (Soderberg et al. 2010a), may be such an event.

Another event, the Type Ic supernova SN 2007gr, is much more controversial.
Soderberg et al. (2010b) proposed that SN 2007gr is an ordinary Type Ic
supernova with an almost constant expansion speed ($v\propto t^{-0.1}$). On
the contrary, Paragi et al.’s (2010) 5 GHz Radio observations have revealed a
relativistic jet in SN 2007gr. While the opening angle of the jet is similar
to that of a typical GRB jet, its Lorentz factor seems to be far smaller than
a normal GRB outflow. In view of the non-relativistic expansion of the
photosphere of SN 2007gr (Valenti et al. 2008; Hunter et al. 2009), here, we
propose an off-axis relativistic jet model for SN 2007gr with typical
parameters of circumstellar medium (CSM) of a Wolf-Rayet star. Most of the
energy in the explosion is contained in non-relativistic ejecta which produces
the supernova, while only a small fraction of expansion energy is contained in
the relativistic jet. We first describe our dynamical model in Section 2. In
section 3 we describe the parameters used in our modeling. The model results
are shown in Section 4. Our conclusions and discussion are presented in
Section 5.

##  2 Model

On 2007 August 15.51 UT, SN 2007gr was discovered by the Katzmann Automatic
Imaging Telescope in NGC 1058 (Madison & Li 2007), a bright spiral galaxy
belonging to a group of galaxies. The distance of one member of this group,
NGC 935, has been derived as about 9.3 Mpc (Silbermann et al. 1996). In our
study, we adopt this value for the distance of SN 2007gr. The explosion date
was suggested as 2007 August $13\pm 2$ (Soderberg et al. 2010b). Radio
observations with the Very Large Array (Soderberg et al. 2010b) and European
VLBI network (EVN) (Paragi et al. 2010) revealed a radio source in the place
of the supernova. SN 2007gr was classified as Type Ic supernova according to
its spectra (Crockett et al. 2008; Valenti et al. 2008).

###  2.1 Relativistic jet

Paragi et al. (2010) reported the evidence for the existence of a relativistic
jet in 2007gr (with expansion speed $v>0.6\rm c$) based on their radio
observations, while Valenti et al. (2008) and Hunter et al. (2009) measured a
non-relativistic velocity (from $\sim 11000~{}\rm km/s$ at 1 week after
explosion to $\sim 4800~{}\rm km/s$ at 50 days after explosion) for the
photospheric expansion of SN 2007gr. These observations suggest that the
average speed of the optical ejecta (i.e., optical expansion velocity) is non-
relativistic while the radio ejecta is much faster and relativistic. This
scenario is similar to that of Soderberg et al. (2005,2010b), but note that
the velocity of the radio ejecta here is relativistic. So, we suggest that SN
2007gr may phenomenologically contain two components: a non-relativistic
component and a relativistic component. The non-relativistic component should
contain most of the explosion energy and account for the photospheric
expansion and the optical emission. The relativistic component should be a jet
that contains only a small fraction of the explosion energy and account for
the radio emission. It is more likely central engine driven. Its behavior
should be similar to a GRB jet, but note that its initial Lorentz factor is
significantly smaller. In this paper, we will simulate the evolution of the
relativistic jet numerically and compare the results with observations.

In our framework, the optical emission of supernova should mainly come from
the decay process of $\rm{}^{56}Ni$ synthesized in the SN explosion. On the
contrary, radio and X-ray emission of supernova is explained as synchrotron
radiation from relativistic electrons which are accelerated by the shock
produced in the collision between the jet and circumstellar medium (Chevalier
1998; Soderberg et al. 2005). The shock process is very similar to the
external shock process of GRBs that gives birth to GRB afterglows.

###  2.2 Jet dynamics

The dynamical evolution of a relativistic jet that collides with surrounding
medium can be conveniently described by the equations proposed by Huang et al.
(1999, 2000). Their method can be widely used in both ultra-relativistic and
non-relativistic phases. In this study, we will adopt Huang et al.’s equations
to simulate the evolution of the external shock. The evolution of the bulk
Lorentz of the jet ($\gamma$), the swept-up mass of medium ($m$), the radius
of the shock ($R$), and the half-opening angle of the jet ($\theta$) are
described by the following equations,

$\frac{d\gamma}{dm}=-\frac{\gamma^{2}-1}{M_{\rm ej}+\varepsilon
m+2(1-\varepsilon)\gamma m},$ (1)

$\frac{dm}{dR}=2\pi R^{2}(1-\cos\theta)nm_{\rm p},$ (2)

$\frac{dR}{dt}=\beta c\gamma(\gamma+\sqrt{\gamma^{2}-1}),$ (3)

$\frac{d\theta}{dt}=\frac{c_{\rm s}(\gamma+\sqrt{\gamma^{2}-1})}{R},$ (4)

where $M_{\rm ej}$ is the initial mass of the ejecta, $m_{p}$ is the mass of
the proton, and $\beta=\sqrt{\gamma^{2}-1}/\gamma$. The radiative efficiency
($\varepsilon$) is assumed as zero because the ejecta becomes adiabatic a few
hours after the burst.

In the above equations, the velocity of the lateral expansion has been assumed
to be the sound speed $c_{\rm s}$, which can be calculated from (Dai et al.
1999)

$c_{\rm
s}^{2}=\frac{\hat{\gamma}(\hat{\gamma}-1)(\gamma-1)c^{2}}{1+\hat{\gamma}(\gamma-1)},$
(5)

where $\hat{\gamma}\approx(4\gamma+1)/(3\gamma)$ is the adiabatic index. The
number density of the circumstellar medium ($n$) is inversely proportional to
the square of the shock radius, i.e.

$n=\frac{\dot{M}}{4\pi m_{\rm p}v_{\rm w}R^{2}},$ (6)

where $\dot{M}$ is mass loss rate of the circumstellar wind, and $v_{\rm w}$
is the wind speed.

###  2.3 Synchrotron radiation process

As usual, in the comoving frame, the shock-accelerated electrons are assumed
to follow a power-law distribution according to their Lorentz factors
($\gamma_{e}$)

$\frac{dN_{\rm e}^{\prime}}{d\gamma_{\rm e}}=(\gamma_{\rm e}-1)^{-p},$ (7)

where $p$ is the power-law index. Note that in the bracket, 1 is subtracted
from $\gamma_{e}$ to account for the non-relativistic phase (Huang & Cheng
2003). For a single electron with a Lorentz factor of $\gamma_{e}$, the
synchrotron radiation power at frequency $\nu^{\prime}$ is given by

$P(\nu^{\prime},\gamma_{\rm e})=\frac{\sqrt{3}e^{3}B^{\prime}}{m_{\rm
e}c^{2}}F(\frac{\nu^{\prime}}{\nu_{\rm c}^{\prime}}),$ (8)

where $e$ is the electron charge, $B^{\prime}$ is the magnetic intensity,
$m_{\rm e}$ is the mass of electron, $\nu_{c}^{\prime}=3\gamma_{\rm
e}^{2}eB^{\prime}/(4\pi m_{\rm e}c)$ and the function $F(x)$ is defined as

$F(x)=x\int_{x}^{+\infty}K_{5/3}(k)dk,$ (9)

with $K_{5/3}(k)$ being the Bessel function.

The magnetic energy density is assumed to be a fraction $\epsilon_{B}^{2}$ of
the total thermal energy density, i.e. (Dai et al. 1999),

$\frac{B^{\prime 2}}{8\pi}=\epsilon_{\rm
B}^{2}\frac{\hat{\gamma}\gamma+1}{\hat{\gamma}-1}(\gamma-1)nm_{\rm p}c^{2}.$
(10)

Therefore, the total synchrotron radiation power from all the shock
accelerated electrons is

$P(\nu^{\prime})=\int_{\gamma_{\rm e,min}}^{\gamma_{\rm e,max}}\frac{dN_{\rm
e}^{\prime}}{d\gamma_{\rm e}}P(\nu^{\prime},\gamma_{\rm e})d\gamma_{\rm e},$
(11)

where $\gamma_{\rm e,max}=10^{8}(B^{\prime}/1~{}\rm G)^{1/2}$ is the maximum
Lorentz factor of elections, $\gamma_{\rm e,min}=\epsilon_{\rm
e}(\gamma-1)m_{\rm p}(p-2)/[m_{\rm e}(p-1)]+1$ is the minimum Lorentz factor
of elections, and $\epsilon_{e}$ is electron energy faction. Then we can
obtain the observed flux density at frequency $\nu$,

$F_{\nu}=\frac{1}{\gamma^{3}(1-\beta\cos\Theta)^{3}}\frac{1}{4\pi
D_{L}^{2}}P^{\prime}[\gamma(1-\beta\cos\Theta)\nu],$ (12)

where $D_{L}$ is the luminosity distance, and $\Theta$ is the angle between
the line of sight and the velocity of emitting material.

The synchrotron self absorption effect (Rybicki & Lightman 1979) need to be
considered in calculating the radio flux (Chevalier 1998; Kong et al. 2009).
Self-absorption reduces the synchrotron radiation flux by a factor of
$(1-e^{\tau_{\nu}})/\tau_{\nu}$, where $\tau_{\nu}$ is the optical depth. The
self-absorption coefficient is given by

$k(\nu^{\prime})=\frac{p+2}{8\pi m_{\rm e}\nu^{\prime 2}}\int_{\gamma_{\rm
e,min}}^{\gamma_{\rm e,max}}\frac{dN_{\rm e}^{\prime}}{d\gamma_{\rm
e}}\frac{1}{\gamma_{\rm e}}P(\nu^{\prime},\gamma_{\rm e})d\gamma_{\rm e}.$
(13)

When calculating the observed flux, we integrate the emission over the whole
equal arrival time surface (Waxman 1997; Sari 1997; Panaitescu & Mészáros
1998) determined by

$t=\int\frac{1-\beta\Theta}{\beta c}dR\equiv{\rm const.}$ (14)

##  3 Parameters

Before fitting the observations of SN 2007gr numerically, we first give a
rough estimation of several fundamental parameters of the jet based on
theoretical analysis of some observed facts.

(i) The initial bulk Lorentz factor $\gamma_{0}$. Paragi et al. (2010) have
derived a conservative lower limit of 1.7 mas for the angular diameter from
their observations made with the EVN and the Green Bank Telescope during 2007
November 5 — 6 ($\sim 85$ days after the supernova explosion). This will
correspond to an average expansion speed faster than 0.6 times the speed of
light ($\beta>0.6$). In view of the deceleration of the ejecta, we assume that
the initial bulk Lorentz factor is 3, i.e., $\beta_{0}\approx 0.94$.

(ii) The initial half-opening angle $\theta_{0}$. The typical half-opening
angle of GRB jets is about $0.1$ rad (Zhang 2007; Gao & Dai 2010), while
supernova outflows are much more isotropic. Since the Lorentz factor of the
jet involved in SN 2007gr is much smaller than a typical GRB jet, its opening
angle should be correspondingly much larger. Here we assume that the initial
half-opening angle of the jet in our model is $\theta_{0}=0.6$ rad ($\sim
34.4^{\circ}$).

(iii) The observing angle $\theta_{\rm obs}$, i.e. the angle between the line
of sight and the symmetry axis of the jet. The ten hours of EVN observations
of SN 2007gr during 2007 November 5 — 6 restored an elliptical beam with the
size of $15.26\times 6.85$ mas (Paragi et al. 2010). In view of the none
detection of prompt emission in the early stage of SN 2007gr, we suggest that
the jet is not pointed toward us, but off-axis. In this paper, the inclination
angle between the jet axis and our line of sight is taken as $\theta_{\rm
obs}=1.1~{}\rm rad$ ($\sim 63.3^{\circ}$). Such an observing angle, together
with the smaller initial half-opening angle of $\theta_{0}=0.6~{}\rm rad$
($\sim 34.4^{\circ}$), means that our line of sight is completely outside the
jet boundry, i.e. we are observing the mildly relativistic jet off-axisly.

(iv) The initial kinetic energy of the shock $E$. The total kinetic energy and
ejected mass of the SN ejecta is $\sim 1.5-3\times 10^{51}$ erg and $\sim
1.5-3.5~{}M_{\odot}$ respectively (Valenti et al. 2008; Hunter et al. 2009).
The total internal energy for the observed radio emission is $\sim
0.7-4.5\times 10^{46}$ erg (Soderberg et al. 2010b; Paragi et al. 2010).
Because most of the energy in the explosion is contained in non-relativistic
ejecta and only very little energy is contained in the relativistic jet, we
assume that the initial kinetic energy of the relativistic shock wave is
$E=1.1\times 10^{48}$ erg.

(v) The number density of the circumstellar medium ($n$). SN 2007gr is
classified as Type Ic SN according to its spectra (Crockett et al. 2008;
Valenti et al. 2008). The most favored progenitors for Type Ic supernovae are
Wolf-Rayet stars (Maeder & Lequeux 1982; Begelman & Sarazin 1986; Woosley &
Bloom 2006). For Wolf-Rayet stars, the typical mass loss rate is $\dot{M}\sim
0.6-17.1\times 10^{-5}~{}M_{\odot}~{}\rm yr^{-1}$, and the wind speed is
typically $v_{\rm w}\sim 0.7-5.5\times 10^{3}~{}\rm km~{}s^{-1}$ (Eenens &
Williams 1994; Cappa et al. 2004). In our calculations, we assume the
following parameters for the progenitor of SN 2007gr, i.e., $\dot{M}=1.0\times
10^{-5}~{}M_{\odot}~{}\rm yr^{-1}$ and $v_{\rm w}\sim 1.5\times 10^{3}~{}\rm
km~{}s^{-1}$. Then the number density of CSM can be calculated from Eq. (6).
The mass loss rate in our modeling is more than 10 times higher than the ones
in previous studies (Paragi et al. (2010); Soderberg et al. (2010b)).

(vi) The electron energy faction $\epsilon_{\rm e}$. In our model, we assume
an evolving electron energy faction as $\epsilon_{\rm e}=\epsilon_{\rm
e,0}\times(R/R_{0})^{\alpha}$. The parameter $\epsilon_{\rm e}$ is evolving
with time. It is slightly different from a normal GRB model. In our
calculations, we assume $\epsilon_{\rm e,0}=0.1$, $R_{0}=3.0\times 10^{16}$
cm, $\alpha=5/4$, i.e.,

$\epsilon_{\rm e}=0.1\times(\frac{R}{3.0\times 10^{16}~{}\rm cm})^{5/4}.$ (15)

Note that our $\epsilon_{\rm e}$ is smaller than unity during all the
observing time.

For other parameters such as the magnetic energy fraction ($\epsilon_{\rm B}$)
and the power-law index of the energy distribution function of electrons
($p$), we take $\epsilon_{\rm B}=0.1$ and $p=3.3$, respectively, which are
similar to those adopted by Soderberg et al. (2010b).

##  4 Results

Using our off-axis relativistic jet model and the parameters above, we have
numerically calculated the multi-band emission of the jet in SN 2007gr. Here
we compare the theoretical light curves with observations.

Fig. 1 shows the light curves of radio emission. From this figure, we see that
the theoretical curves at 1.4 GHz (Fig. 1a), 4.9 GHz (Fig. 1b), and 8.5 GHz
(Fig. 1c) can fit the observational data well (note the last data point in
Fig. 1c is an upper limit). The curve at 22.5 GHz (Fig. 1d) is also consistent
with the upper limit of the observation.

Figure 1: Observed radio light curves of SN 2007gr and our best fit by using
the off-axis relativistic jet model. Panels (a), (b), (c), and (d) are light
curves at 1.4, 4.9, 8.5 and 22.5 GHz respectively. Observed data points are
taken from Soderberg et al. (2010b).

Fig. 2 illustrates the theoretical optical and X-ray light curves based on our
model. Our X-ray light curve is consistent with the upper limit of the
$Chandra$ observation (Soderberg 2007). The predicted optical emission is
significantly lower than the observed flux. This is a reasonable result. We
believe that the observed optical emission should mainly come from the decay
process of $\rm{}^{56}Ni$ synthesized in the bulk SN ejecta (Arnett 1982;
Sutherland & Wheeler 1984), as that typically happens in usual SN explosions.

Figure 2: Light curves of SN 2007gr in optical and X-ray bands. The solid and
dashed curves correspond to optical and X-ray emission from our off-axis
relatisvistic jet, respectively. The circular points are the observed R-band
data taken from Valenti et al. (2008). The triangle point is the upper limit
observed by $Chandra~{}X-ray~{}Observatory$ (Soderberg 2007).

Figs. 1 and 2 are our best fits to the observations by adopting optimal values
for the parameters involved. To get these optimal parameters, we actually have
tried many times. In Figs. 3 and 4, we go further to give some examples of the
radio light curves to illustrate the effects of various parameters. In these
figures, the observed 1.4GHz data points are taken from Soderberg et al.
(2010b). The solid curve corresponds to our best fit by using the parameters
described in Section 3. The dashed and dotted curves are drawn with only one
parameter altered. From these figures, we see that the theoretical radio
fluxes depend sensitively on the parameters of $\theta$, $\dot{M}$, and
$v_{\rm w}$. The flux of a more isotropic outflow is higher than that of a
narrower outflow. A fast wind with a low mass loss rate tends to increase the
emission. The parameters $\gamma,~{}E,~{}\epsilon_{\rm B}$, and $\epsilon_{\rm
e}$ affect not only the intensity, but also the shape and peak time of the
light curve. It means that the radiation spectra depends on these parameters
sensitively.

As is shown in Fig. 4b, the light curve is flattened by a positive $\alpha$.
If $\epsilon_{\rm e}$ is constant as normally assumed in GRB modeling, i.e.,
$\alpha=0$, then the flux decreases much faster than the positive $\alpha$
scenario. A positive $\alpha$ is necessary in our fitting, which indicates
that $\epsilon_{\rm e}$ is increasing with time, i.e., more and more kinetic
energy is gradually transformed to internal energy.

Figure 3: Effects of the parameters $\gamma$, $\theta$, $E$, and
$\epsilon_{\rm B}$ on the 1.4 GHz radio light curve. Observed data points of
SN 2007gr are again taken from Soderberg et al. (2010b). The solid curves are
our best fit, while the dashed and dotted curves are light curves with only
one parameter altered, as marked in the panels.  Figure 4: Effects of the
parameters $\epsilon_{\rm e}$, $\dot{M}$, and $v_{\rm w}$ on the 1.4 GHz radio
light curve. Observed data points and the solid curves are the same as in Fig.
3. The dashed and dotted curves are light curves with only one parameter
altered, as marked in the panels.

##  5 Discussion and Conclusions

In this paper, we propose an off-axis relativistic jet model to explain the
multiband emission of the Type Ic supernova SN 2007gr. From the observations
of SN 2007gr by EVN (Paragi et al. 2010), we adopted a model where the outflow
contains a jet component and the line of sight is off-axis. The observing
angle is about $\theta_{\rm obs}=1.1$ rad. We suggest that the radio emission
come from synchrotron radiation of relativistic electrons accelerated by the
shock produced by the collision between the jet and the circumstellar medium.
Our calculations show that a jet with an initial half-opening angle of
$\theta_{0}=0.6$ rad and initial Lorentz factor of $\gamma_{0}=3$ can
reproduce the emission in radio band well. Optical and X-ray emission is also
consistent with the observational restrictions.

The jet engaged in our modeling of SN 2007gr differs significantly from normal
GRB jet. First, the Lorentz factor of normal GRB ejecta is typically several
hundred, while it is only a few here for SN 2007gr. Second, the opening-angle
of our SN jet ($0.6~{}\rm rad$) is much larger than that of a GRB jet ($\sim
0.1~{}\rm rad$). Additionally, we would like to point out that the lateral
expansion speed of the jet is a relatively complicate factor. In our
calculations, we have assumed it to be the comoving sound speed. Although this
is a reasonable assumption, deviation may still be possible and the dynamics
may be affected. Third, the initial energy of our SN jet (isotropically $\sim
10^{48}~{}\rm erg$) is much smaller than that of GRB jet ($\sim 10^{52}~{}\rm
erg$). We believe that the majority of energy is deposited into supernova
component in SN 2007gr explosion. Thus the above characteristics are not too
difficult to understand.

In view of Wolf-Rayet stars as the progenitors of type Ic SNe, we assume the
CSM as the typical wind of a Wolf-Rayet star with mass-loss rate
$\dot{M}=1.0\times 10^{-5}~{}M_{\odot}~{}\rm yr^{-1}$ and $v_{\rm w}=1.5\times
10^{3}~{}\rm km~{}s^{-1}$. A fast wind with low mass loss rate tend to
increase the flux of the radio emission. The mass loss rate in our model is
more than 10 times higher than that in previous studies, but it is more
typical for a Wolf-Rayet star. Note that in our modeling, we have assumed that
the CSM is purely composed of protons and electrons. Actually, other nucleus
such as helium may also appear in the wind of Wolf-Rayet stars. Although their
presence in the CSM should not change the final results significantly, we may
still need to consider the factor in the future, when observations become much
improved.

For the electron energy faction ($\epsilon_{\rm e}$), we usually assume that
it is constant ($\alpha=0$) in GRB model. However, in the current study, a
positive $\alpha$ is necessary in our fitting. It means that $\epsilon_{\rm
e}$ is increasing with time, i.e., the fraction of energy going to internal
energy increases with the deceleration of the shock. Although the assumption
of a varying $\epsilon_{\rm e}$ has also been engaged in modeling some special
GRBs and other transient objects (Rossi & Rees 2003; Ioka et al. 2006; Kong et
al. 2010), note that the underlying physical mechanism that leads to the
variation of $\epsilon_{\rm e}$ is still much uncertain.

The mutiband emission and spectra are quite different between SN 2007gr and
normal GRBs. As for GRBs, we can detect the prompt gamma-ray emission and
multiband afterglows. The peak frequency of prompt emission of GRBs is
typically in gamma-ray band (about several hundred keV). The peak frequency of
GRB afterglows is typically in UV or X-ray band, and the X-ray afterglow can
usually be observed. However, the emission from the jet of SN 2007gr is more
like a failed GRB rather than a normal GRB, since it is only a midly
relativistic outflow. It has been proposed that the prompt emission of a midly
relativistic SN jet come from the photosphere and bright in UV or soft X-ray
band (Xu et al. 2010). The peak frequency of afterglow should then be in near
infrared or radio band. The X-ray emission should be weak and hard to detect.

In the future, more and more relativistic SN ejecta will be detected. These
relativistic ejecta are considered as central engine driven. In fact,
supernova 2009bb was observed as another example of relativistic supernova
ejecta recently (Soderberg et al. 2010b). If the very early emission of
supernovae could be detected in the future, it will be helpful for determining
the speed of the supernova ejecta directly.

We thank the anonymous referee for helpful comments and suggestions. This work
was supported by the National Natural Science Foundation of China (Grant No.
10625313 and 11033002), the National Basic Research Program of China (973
Program, Grant No. 2009CB824800) and the Grant-in-Aid for the ‘Global COE
Bilateral International Exchange Program’ of Japan, Grant-in-Aid for
Scientific Research on Priority Areas No. 19047004 and Scientific Research on
Innovative Areas No. 21105509 by Ministry of Education, Culture, Sports,
Science and Technology (MEXT), Grant-in-Aid for Scientific Research (S) No.
19104006 and Scientific Research (C) No. 21540404 by Japan Society for the
Promotion of Science (JSPS).

## References

  * Arnett (1982) Arnett, W. D. 1982, ApJ, 253, 785
  * Barkat et al. (1974) Barkat, Z., Wheeler, J. C., Buchler, J.-R., & Rakavy, G. 1974, Ap&SS, 29, 267
  * Begelman & Sarazin (1986) Begelman, M. C., & Sarazin, C. L. 1986, ApJ, 302, L59
  * Beger et al. (2003) Berger, E. et al. 2003, Nature, 426, 154
  * Campana et al. (2006) Campana, S. et al. 2006, Nature, 442, 1008
  * Cappa et al. (2004) Cappa, C., Goss, W. M., & van der Hucht, K. A. 2004, AJ, 127, 2885
  * Chevalier (98) Chevalier, R. A. 1998, ApJ, 499, 810
  * Cobb et al. (2004) Cobb, B. E., Bailyn, C. D., van Dokkum, P. G., Buxton, M. M., & Bloom, J. S. 2004, ApJ, 608, L93
  * Cobb et al. (2010) Cobb, B. E., Bloom, J. S., Perley, D. A., Morgan, A. N., Cenko, S. B., & Filippenko, A. V. ApJ, 2010, 718, 150
  * Crockett et al. (2008) Crockett, R. M. et al. 2008, ApJ, 672, L99
  * Dai et al. (1999) Dai, Z. G., Huang, Y. F., & Lu, T. 1999, ApJ, 520, 634
  * Eenens & Williams (1994) Eenens, P. R. J., & Williams, P. M. 1994, MNRAS, 269, 1082
  * Fan et al. (2010) Fan, Y. Z., Zhang, B. B., Xu, D., Liang, E. W., & Zhang, B. 2010, ApJ, 762, 32
  * Filippenko (1997) Filippenko, A. V. 1997, ARA&A, 35, 309.
  * Frail et al. (2001) Frail, D. A. et al. 2001. ApJ, 562, L55
  * Gehrels (2009) Gehrels, N., Ramirez-Ruiz, E., & Fox, D. B. 2009, ARA&A, 47, 567
  * Galama et al. (1998) Galama, T. J. et al. 1998, Nature 395, 670
  * Gao & Dai (2010) Gao, Y., & Dai, Z.G. 2010, RAA, 10, 142
  * Granot & Loeb (2003) Granot, J., & Loeb, A. 2003, ApJ, 593, L81
  * Huang & Cheng (1999) Huang, Y. F., Cheng, K. S. 2003, MNRAS, 341, 263
  * Huang et al. (1999) Huang, Y. F., Dai, Z. G., & Lu, T. 1999, MNRAS, 309, 513
  * Huang et al. (2000) Huang, Y. F., Gou, L. J., Dai, Z. G., & Lu, T. 2000, ApJ, 543, 90
  * Hunter et al. (2009) Hunter, D. J. et al. 2009, A&A, 508, 371
  * Ioka et al. (2006) Ioka, K., Toma, K., Yamazaki, R., & Nakamura, N. 2006, A&A, 458, 7
  * Kong et al. (2009) Kong, S. W., Huang, Y. F., Cheng, K. S., & Lu T. 2009, Science in China Series G, 52, 2047
  * Kong et al. (2010) Kong, S. W., Wong, A. Y. L., Huang, Y. F., & Cheng, K. S. 2010, MNRAS, 402, 409
  * Madison & Li (2007) Madison, D., & Li, W. 2007, CBET, 1034, 1
  * Maeder & Lequeux (1982) Maeder, A., & Lequeux, J. 1982, A&A, 114, 409
  * Mészáros (2002) Mészáros, P. 2002, ARA&A, 40,137
  * Nagataki et al. (2007) Nagataki, S., Takahashi, R., Mizuta, A., Takiwaki, T. 2007, ApJ, 659, 512
  * Nagataki (2009) Nagataki, S. 2009, ApJ, 704, 937
  * Nagataki (2010) Nagataki, S. 2010, PASJ, submitted (arXiv:1010.4964)
  * Paragi et al. (2010) Paragi, Z. et al. 2010, Nature, 463, 516
  * Rossi & Rees (2003) Rossi, E., & Rees, M. J. 2003, MNRAS, 339, 881
  * Rybicki & Lightman (1979) Rybicki, G. B., & Lightman, A. P. 1979, _Radiative Processes in Astrophysics_ , New York: Wiley
  * Silbermann et al. (1996) Silbermann, N. A. et al. 1996, ApJ, 470, 1
  * Soderberg (2007) Soderberg, A. M. 2007, ATel, 1205, 1
  * Soderberg et al. (2005) Soderberg, A. M. et al. 2005, ApJ, 621, 908
  * Soderberg et al. (2006) Soderberg, A. M., Nakar, E., Berger, E., & Kulkarni, S. R. 2006, ApJ, 638, 930
  * Soderberg et al. (2010) Soderberg, A, M, et al. 2010a, Nature, 463, 513
  * Soderberg et al. (2010) Soderberg, A. M., Brunthaler, A., Nakar, E., & Chevalier, R. A. 2010b, ApJ, 725, 922
  * Sutherland & Wheeler (1984) Sutherland, P. G., & Wheeler, J. C. 1984, ApJ, 280, 282
  * Valenti et al. (2008) Valenti, S. et al. 2008, ApJ, 673, L155
  * Wheeler (1993) Wheeler, J. C. 1993, PKAS, 8, 169
  * Wheeler & Levreault (1985) Wheeler, J. C., & Levreault, R. 1985, ApJ, 294, L17
  * Wheeler et al. (1993) Wheeler, J. C., Swartz, D. A., & Harkness, R. P. 1993, PhR, 227, 113
  * Wilson (1971) Wilson, J. R. 1971, ApJ, 163, 209
  * Woosley & Bloom (2006) Woosley, S. E., & Bloom, J. S. 2006, ARA&A, 44, 507
  * Woosley & Janka (2005) Woosley, S. E., & Janka, T. 2005, Nat. Phys., 3,147
  * Xu et al. (2010) Xu, M, Nagataki, S., & Huang, Y. F. 2010, ApJ submitted
  * Zhang (2007) Zhang, B. 2007, ChJAA, 7, 1