# Current rectification and ionic selectivity of $\alpha$-hemolysin: Coarse-
Grained Molecular Dynamics simulations

Delphine Dessaux  Jérôme Mathé  Rosa Ramirez  Nathalie Basdevant
<EMAIL_ADDRESS>Université Paris-‐Saclay, CNRS, Univ Evry,
LAMBE, 91025, Évry‐-Courcouronnes, France

###### Abstract

In order to understand the physical processes of nanopore experiments at the
molecular level, microscopic information from molecular dynamics is greatly
needed. Coarse-grained models are a good alternative to classical all-atom
models since they allow longer simulations and application of lower electric
potentials, closer to the experimental ones. We performed coarse-grained
molecular dynamics of the ionic transport through the $\alpha$-hemolysin
protein nanopore, inserted into a lipid bilayer surrounded by solvent and
ions. For this purpose, we used the MARTINI coarse-grained force field and its
polarizable water solvent (PW). Moreover, the electric potential difference
applied experimentally was mimicked by the application of an electric field to
the system. We present, in this study, the results of 1.5 microsecond long-
molecular dynamics simulations of twelve different systems for which different
charged amino acids were neutralized, each of them in the presence of nine
different electric fields ranging between $\pm 0.04$ V.nm-1 (a total of around
100 simulations). We were able to observe several specific features of this
pore, current asymmetry and anion selectivity, in agreement with previous
studies and experiments, and identified the charged amino acids responsible
for these current behaviors, therefore validating our coarse-grain approach to
study ionic transport through nanopores. We also propose a microscopic
explanation of these ionic current features using ionic density maps.

###### keywords:

Molecular Dynamics, Coarse-Grain, $\alpha$-hemolysin, nanopore, simulations,
ionic current

MD, CG, AA, $\alpha$HL

##  1 Introduction

The nanopore technique is a popular method to study the kinetics and
energetics of biological phenomena. It consists in applying an electric
potential difference to guide a charged polymer through an artificial or
protein nanopore, inserted in a solid or lipid membrane, in the presence of a
salt solution 1, 2, 3, 4, 5. A macromolecule passing through the pore induces
a decrease in the ionic current by partially blocking it, depending on the
nature of the molecule and the pore characteristics: size, conformation, and
structure. The $\alpha$-hemolysin ($\alpha$HL) is a biological protein
nanopore widely employed in nanopores experiments. Inserted in a lipid
bilayer, it has been used to study DNA or RNA translocation 1, 3, DNA
sequencing or biosensing 6, 7, 8, DNA unzipping of individual hairpins 9, 10,
11, and the translocation of proteins 4, 12, 13, 14, polysaccharides 15 or
synthetic polymers 16, 17.

Several experimental 1, 18, 19, 20, theoretical 21 and computational studies
22, 23, 24, 25, 26, 27, 28, 29 on the ionic current through $\alpha$HL
inserted in a lipid bilayer were already performed and they have shown that
this current is not symmetrical for positive and negative applied electric
potentials. This effect is called current rectification. It can be explained
by the asymmetric charge distribution of the amino acids of the pore which
yields an asymmetric electric potential, resulting in an asymmetric crossing
of ions under a positive or negative applied voltage difference. It should be
noted that this effect is specific to each protein pore. Indeed, the
rectification observed for other protein nanopores differs from $\alpha$HL, as
exemplified by Payet et al. for the aerolysin pore 20. Zhou et al. 29 recently
showed with MD simulations of several protein nanopores that the channel
geometry and charge distribution affect ion transport. In addition, several
experiments pointed that cation type, temperature, ionic concentration and pH
can influence the rectification of $\alpha$HL. For example, Piguet et al. 19
showed that the rectification is more important for KCl solution than for
LiCl. Battacharya et al. 25 performed experiments and Molecular Dynamics (MD)
simulations with different cation types and found a greater conductance for K+
than for Na+ and Li+, which they associated to the cations ability to screen
the pore amino acids charges. Besides, Payet et al. 20 conducted experiments
showing that, while the ionic current increases with temperature, the current
asymmetry decreases. Several ionic transport experimental studies 30, 18, 31
with different KCl concentrations have also proven that conductance increases
with concentration whereas asymmetry decreases. As for the pH influence, it
has been shown both by experiments 21 and all-atom MD simulations 28 that the
current asymmetry decreases with low pH. In this study, we will focus on the
current asymmetry and ionic selectivity of $\alpha$HL.

The $\alpha$HL nanopore structure counts seven identical chains composed of
293 amino acids, and it involves a hydrophobic stem part (mostly composed of
beta-sheets), which is embedded into the lipid membrane, and a cap part,
facing the cis side of the membrane (see Fig. 1). The trans part of $\alpha$HL
includes three charged amino acids on each of the seven chains: two aspartic
acids (D127 and D128) and one lysine (K131), resulting in a total charge of
-7e (see Figure 1). The lateral chains of D127 and K131 point at the inner
channel, while D128 are directed towards the membrane. Studies involving
Poisson-Nernst-Planck (PNP) theory 22, as well as all-atom MD simulations 25,
revealed that the neutralization of the charged amino acids at the trans part
induces the disappearance of the current rectification. Moreover, Merzlyak et
al. 18 performed mutations of some charged amino acids in the pore stem which
revealed that the residue charges at the trans part of the pore affect the
rectification.

Furthermore, the $\alpha$HL pore is slightly anion selective. This means that,
for all applied voltages through the pore, more anions are crossing over than
cations. The charged amino acids inside the protein stem also play a role in
this effect. In addition to the three charged amino acids at the trans part,
there are two charged amino acids at the top of the stem, i.e. at the pore
constriction: one glutamic acid (E111) and one lysine (K147) on each of the
seven chains. Neutralizing these two amino acids affects the anion selectivity
of the pore. Using PNP 22, these two charged residues were indeed identified
to influence the selectivity. Besides, experiments and all-atom MD 25 showed
that this anion selectivity was more important with K+ cations than with Na+
or Li+. However, both theoretical studies pointed out that neutralization of
the charged residues at the _trans_ part increased weakly the selectivity.
Moreover, mutation experiments 32, 33 revealed that the mutation of E111 and
K147 to neutral residues results in cation selectivity. Figure 1 represents
the $\alpha$HL structure (PDB: 7AHL) and the location of the charged amino
acids in the stem.

Figure 1: Representation of $\alpha$-hemolysin in cartoon mode. The five
charged amino acids in the stem are highlighted in different colors. At the
pore constriction (E111 and K147) and at the trans side (D127, D128 and K131).
On the right, view from the trans side.

Here, we choose to study these two effects, rectification and selectivity,
performing Coarse-Grained (CG) molecular dynamics simulations, using the
MARTINI force field. We previously showed that this approach is adapted to
study the ionic current through $\alpha$HL inserted into a DPPC bilayer at 1 M
KCl concentration 34. Indeed, CG-MD simulations allow to simulate the system
for longer timescales than All-Atom (AA) simulations and with voltages closer
to experimental ones. Though, we showed that there is a scaling factor for
ionic conductances: the simulated currents are lower than experimental ones.
In this study, we aim to identify the role of the charged amino acids in
current asymmetry and anion selectivity. Therefore, we calculated ionic
currents in the presence of several external electric fields, when specific
charged amino acids inside the pore stem were neutralized, performing around
100 different simulations of 1.5 microsecond. We decomposed the ionic current
into its anion and cation contribution, in order to evaluate the pore
selectivity. We tested many more different combinations of amino acid
neutralization compared to previous AA-MD simulations. We were thus able to
identify the amino acids responsible for the rectification, at the trans side
of the pore, and the ones responsible for the anion selectivity, at the
constriction part of the stem. We explained these effects, using ionic density
maps, comparing the anions and cations behavior in the vicinity of the pore
surface.

##  2 Results and discussion

All our CG-MD simulations are performed following the same protocol (see
Methods) as in our previous study 34: we use the MARTINI CG force field for
the $\alpha$HL protein inserted into a DPPC bilayer, surrounded by solvent (PW
polarizable CG water molecules) and NaCl ions at a 1 M concentration. In
addition, ElNeDyn elastic network 35 was applied on the protein structure. We
present the simulations of several systems for which we neutralized some of
the charged amino acids on the seven chains of the $\alpha$HL protein. We
refer to these simulations using the following numbering: e.g. D127N refers to
the simulations for which the aspartic acid 127 of each chain was neutralized.
Moreover, for each system, in order to mimic the voltage difference applied
experimentally, several electric fields were applied perpendicularly to the
membrane plan, ranging from $-0.04$ V/nm to $+0.04$ V/nm. It was shown that it
is equivalent to applying a potential difference $\Delta V=E_{z}.L_{z}$ 36.
Several previous MD studies on ionic transport through nanopores confirmed
this relation 37, 23, 38, 39.

From all these simulations, we generate the $IV$ curves of each pore and we
calculated their ionic conductances $G^{+}$ and $G^{-}$ using two linear
regressions, one for positive voltages and one for negative voltages. With the
same method, we compute the separate conductances: $G_{\text{Cl}}$ for anions
and $G_{\text{Na}}$ for cations. For all our CG-MD simulations, we evaluate
the current asymmetry using the conductance rectification ratio $r$ between
positive and negative potentials: $r=\frac{G^{+}}{G^{-}}$ (see Methods for
details). We also determine the anion asymmetry ratio:
$r_{\text{Cl}}=\frac{G_{\text{Cl}}^{+}}{G_{\text{Cl}}^{-}}$, in order to
identify whether the rectification effect comes from the anionic or cationic
current. The ionic selectivity is defined by the selectivity ratio $s$ between
the conductances of cations and anions for a positive electric potential:
$s=\frac{{G_{\text{Na}}}^{+}}{{G_{\text{Cl}}}^{+}}$. An $s$ ratio equal to one
stands for a non selective pore, if it is greater than one, the pore is cation
selective, and if it is lower than one, the pore is anion selective. Table 1
compiles rectification and selectivity ratios for all simulations, as well as
the average conductances for positive and negative electric fields. We also
draw the average ionic density maps for all simulations, in order to
understand how charge distribution around the pore affects the ionic flow.
This map is derived from the charged particle density averaged over the 1 to
1.5 $\mu$s part of the simulations, centered on the constriction part (see
Methods).

###  2.1 Wild-type and neutralized $\alpha$-hemolysins

We first compare our results for the native heptameric protein (WT, wild-type)
to a system with a totally neutralized pore, in order to validate our
protocol. The WT pore carries a global charge of $+7$e, corresponding to the
state of the protein in physiological conditions, whereas for the totally
neutral pore, all protein partial charges were set to zero.

Without any restraint on the neutralized $\alpha$HL structure, we find almost
no current passing through the pore. We observe a closure of the channel at
the trans extremity, particularly for low electric fields. We explain this
effect by the lack of elastic network in this region. Since the elastic
network is applied within a cut-off of 0.9 nm, the presence of loops at the
trans part of the channel is responsible for a gap between the different
chains, which does not allow the addition of many springs. Therefore, this
part of the structure is less stabilized by the elastic network. For the
native protein, coulombic repulsion between charged amino acids helps to
maintain the channel open. On the contrary, for the neutral pore, the trans
extremity is unstable and collapses. We observe the same closure of the trans
extremity in other CG-MD simulations with a modified $\alpha$HL at the trans
side (especially D127N and D127N-D128N-K131N ; data not shown). To avoid this
closure effect, for all systems with a pore modified at the trans extremity,
position restraints are applied on the protein CG backbone beads to their
initial positions. Pore simulations using these position restraints are called
thereafter ”constrained pores”, or distinguished by an asterisk (*).

The current asymmetry for the WT pore was previously presented and briefly
discussed 34. It should be noted that, as we have also shown in our last
paper, the CG currents are almost 15 times lower compared to experiments,
which is due to the CG-particles size. Using the computed $IV$ curve of the WT
pore (Figure 2A), we deduce a current asymmetry ratio of $r=1.35$ for the
native pore at 1 M NaCl concentration, 320 K temperature and pH 7 (see Table
1). This value can be interpreted as a 35 $\%$ asymmetry, which is consistent
with previous AA-MD studies 23, 28 and PNP and Brownian dynamics simulations
22. Experiments under the same conditions also show a 30% asymmetry 20. In
Figure 2A, we see that the ionic current through the pore is mostly due to the
chloride anions. However, the chloride conductances for positive and negative
potentials are almost identical, unlike for the sodium cations, as confirmed
by the chloride rectification factor $r_{\text{Cl}}$ in Table 1
($r_{\text{Cl}}=0.91$). Indeed, as shown on the $IV$ curve, the cation
conductance for negative electric fields (${G^{-}_{\text{Na}}}$) is almost
zero. Thereby, the current rectification effect measured for $\alpha$HL comes
from the cationic current through the pore. Moreover, the selectivity ratio
calculated from anion and cation conductances is $s=0.60$, showing that the WT
$\alpha$HL is highly anion selective, as it was previously observed by both
experiments and simulations 30, 40, 32, 22, 18, 33, 41, 25, 42, 19, 43, 28.

(A)0cm-3.9cm (B)0cm-3.9cm

Figure 2: $IV$ curves for the WT pore (A) and neutral constrained $\alpha$HL
(B). We show in green the lines resulting from linear regression either on the
positive or the negative potentials.

We also studied the ionic current behavior for a neutralized constrained pore.
Figure 2B represents the corresponding $IV$ curve. We observe the absence of
current rectification, highlighted by its current asymmetry ratio $r=1.03$,
and a low anion selectivity ($s=0.85$). Indeed, contrary to the WT pore, a
cationic current is measured for negative electric fields. Furthermore,
cationic and anionic currents are almost equal, regardless of the sign of the
electric field. We can also notice that global current intensities are higher
compared to the results for the WT pore with identical applied voltage (see
conductance values in Table 1), which is confirmed by a scaling factor from
experimental currents of only 8 instead of 15 for the WT pore. This is
explained by a greater cationic current for the neutralized pore. These
results are in accordance with the previous theoretical study of Noskov et al.
who obtained a current increase for systems with an $\alpha$-hemolysin
neutralized at the top of the stem region 22.

(A)0.2cm-2.0cm (B)0.2cm-2.5cm

Figure 3: Ionic density maps for the WT $\alpha$-hemolysin, in the presence of
a positive electric field $E_{z}=+0.03$ V/nm (A) and negative electric field
$E_{z}=-0.03$ V/nm (B).

In addition, we plotted the average ionic density map for the WT pore, in
order to locate the regions where the protein pore interacts strongly with
ions. The average density is calculated on the last 500 ns of each trajectory
(see Methods). Figure 3 shows the average density maps for positive (A) and
negative (B) electric fields ($E_{z}=\pm 0.03$ V/nm) of the WT pore. Ions
accumulate at both extremities of the stem, for both positive and negative
potentials. Sodium ions (in blue) are located at both sides, while chloride
ions (in red) are mostly located just above the constriction part. Only a few
ions are found inside the stem in the presence of a negative electric field.
However, with a positive electric field, there are several interaction sites,
mostly with sodium ions. This effect is consistent with the calculation of the
total number of ions in the stem, showing many more ions when a positive
electric field is applied (see Supplemental Information). Furthermore, these
density maps point out that the difference between ionic densities with
positive and negative electric fields mostly comes from the sodium ions. This
result is in agreement with the fact that current rectification is due to the
cationic current.

###  2.2 $\alpha$-hemolysin with modified trans side

It has previously been shown that the charged amino acids at the trans side of
$\alpha$HL are directly involved in the current rectification phenomenon:
Noskov et al. 22 using 3D-PNP and Brownian Dynamics, Bhattacharya et al. 25
using AA-MD simulations, Misakian and Kasianowicz performing pH experiments 21
and Merzlyak et al. 18 conducting mutagenesis experiments. We therefore
computed the ionic current behavior of _trans_ amino acids modified pores.

We performed CG-MD simulations on seven $\alpha$HL with a modified trans part
and, as explained in the previous part of this article, position restraints
are applied on the backbone of the protein for these simulations to prevent
the pore closure.

The results for all our simulations with a modified trans extremity are
compiled in Table 1. We first focus on the charge influence of the aspartic
acids D127 and D128 by neutralizing them separately, resulting in an
$\alpha$HL with a neutral trans part. We simulated thereby two systems: one
with D127 protonated (D127N*) and one with D128 protonated (D128N*). Figure 4A
shows the $IV$ curve obtained for D128N* system. No current rectification is
observed for this pore ($r=1.0$), as well as a strong anion selectivity with a
ratio $s=0.48$. Similarly to the WT pore, the cationic current is asymmetric
and almost null for negative electric fields, but the anionic current shows a
strong reverse asymmetry. The combination of these two currents results in no
rectification. On the other hand, as shown on Table 1, for D127N*, the current
asymmetry is reversed compared to the WT pore, with a ratio $r=0.60$, but the
anion selectivity is still present ($s=0.40$). Our results highlight a
different role for these two amino acids, and are consistent with the studies
of Noskov et al. 22 and Bhattacharya et al. 25 who found no rectification but
a greater anion selectivity when the D128 were neutralized. Previous mutation
experiments 18 showed only a low decrease in rectification when D127 are
mutated into cysteines, and a reverse asymmetry when D128 were mutated into
cysteines. The fact that cysteines are smaller and that their polarity differs
from neutral aspartic acids can explain the discrepancies between these
experiments and our simulations using neutralized aspartic acids.

We also neutralized all trans charged residues (D127, D128 and K131 on each of
the seven monomers: D127N-D128N-K131N*). Table 1 shows that there is almost no
current rectification ($r=1.04$) for this pore, but a huge anion selectivity
with a ratio $s=0.29$. These results are consistent with previous theoretical
studies 22 and AA-MD 25, and this tends to confirm the validity of our
protocol. Our results also suggest that the rectification effect and the anion
selectivity are not totally independent phenomena.

(A)0cm-3.9cm (B)0cm-3.9cm

Figure 4: $IV$ curves for the constrained $\alpha$HL with a modified _trans_
extremity: D128N* (A) and D127N-K131N* (B).

We simulated some additional trans side modified pores: both D127 and D128
protonated (D127N-D128N*, with a trans charge of $+7$e), only K131
deprotonated (K131N*, trans charge of $-14$e), or K131N in combination with
D127N or D128N, therefore with the same trans charge than the WT-pore. We
observe different behavior regarding the current rectification, while all
these pores remained highly anion selective with ratios between $s=0.32$ and
$s=0.49$, as presented in Table 1. Regarding rectification, as shown for
D127N-K131N* in Figure 4B and D128N-K131N* in Table 1, the absence of charge
on K131 in combination with the neutralization of one of the two aspartic
acids results in a reverse asymmetry (respectively $r=0.60$ and $r=0.61$),
although the total trans charge is the same as the WT $\alpha$HL. For these
two systems, this effect comes from a greater anionic current for negative
electric fields compared to WT. For D127N-D128N* and for K131N*, a loss of
current asymmetry was obtained ($r=1.08$ and $r=0.96$, respectively). These
results suggest that this is not the total charge of the trans part of the
stem which is responsible for rectification, but the delicate balance between
the charges of the three charged amino acids which can influence the ionic
distribution in the stem, as will be discussed thereafter.

(A)0.2cm-2.0cm (B)0.2cm-2.0cm (C)0.2cm-2.5cm

Figure 5: Ionic density maps for the modified $\alpha$-hemolysins D128N* (A),
D127N* (B) and D127N-K131N (C) in the presence of a negative electric field
$E_{z}=-0.03$ V/nm.

Figure 5 represents average ionic density maps, in the presence of a negative
electric field of $E_{z}=-0.03$ V/nm (differences in the density maps are more
noticeable with negative fields), for the D128N* (A), D127N* (B) and
D127N-K131N* (C) $\alpha$HL pores. D128N* shows no rectification, while D127N*
and D127N-K131N* present an inverse rectification. Compared to the WT pore, we
observe a higher negative ionic density inside the stem for these three
systems in the presence of a negative electric field. This is consistent with
a null cationic current at negative voltages for all this pores compared to
the already weak cationic current of the WT.

All these density maps are very similar at the constriction part. However,
D128N* (Fig. 5A) and D127N* (Fig. 5B) display very different characteristics
of ionic densities at the trans part, despite sharing the same charge. Indeed,
D127 residues point towards the inner pore and are very close to K131, while
the D128 amino acids point towards the outer pore (see Fig. 1). D127 and K131
tends to form salt bridges, and the neutralization of one of these residues
breaks this bond, therefore affecting the current asymmetry and the ionic
density.

For the D128N* pore, we notice on Fig. 5A two peaks of positive and negative
ionic densities of similar intensity at the trans part. Since the total charge
for this part of the protein is zero, the resulting global charge of the
system at the trans part is almost null. This is consistent with the lack of
current rectification of this pore.

For the D127N* system, we observe on Fig. 5B a highly negative spot at the
trans side of the stem, which is not present for the other pores, including
the WT. The accumulation of anions at this location compensates the positive
charges of the K131 residues. The bottom global charge is therefore negative.

The D127N-K131N* pore on Fig. 5C, bearing the same charge as the WT
$\alpha$HL, does not reveal the positive density peak at the trans part
observed for WT pore (see Figure 3). Thus, the -7$e$ negative charge at the
bottom of the pore is not compensated by the presence of cations, resulting in
a bottom global charge similar to D127N*. The presence of this negative global
charge at the trans part causes the reverse rectification phenomenon. We
should notice that, whereas D127N* and D127N-K131N* systems both exhibit the
same reverse rectification, the currents are twice as low for D127N*. We can
explain this phenomenon by the absence of ion accumulation at this side of the
D127N-K131N* pore which induces a higher conductance compared to D127N*.

The density maps of the other $\alpha$HL pores can be found in the
Supplemental Information. The study of these ionic distributions confirms the
role of the trans part charges on current rectification, and provides an
explanation for the different rectification behaviors of the pores despite
similar total charges.

Our analyses confirm that the three charged amino acids of the $\alpha$HL
trans are responsible for the rectification but play no direct role in the
anion selectivity. We also identified that D127 and D128 play different roles
because of their different positioning. As neutralizing some amino acids at
the trans part (D128N*, D127N-D128N-K131N*, K131N* and D127N-D128N*) cancels
the rectification effect, the neutralization of certain combinations of amino
acids can reverse the current asymmetry (D127N*, D127N-K131N* and
D128N-K131N*). These behaviors can be explained by the distribution of ions
within the channel, apparent in our ionic density maps. All these systems keep
an anion selectivity equivalent or greater than the WT pore.

###  2.3 $\alpha$-hemolysin with modified constriction

In this section, we focus on the charged residues located in the $\alpha$HL
constriction, at the top of stem: E111 and K147 (see Fig. 1). The constriction
is the narrower part of the pore, where the channel diameter reduces from 4.6
to 1.4 nm. It carries a null total charge for the WT pore and, since the two
charged amino acids form a salt bridge, a modification of the charges of these
residues should influence the ionic flow. A theoretical study from Noskov and
Roux 22 showed that the charges at the constriction influence anion
selectivity and rectification. Moreover, several experimental mutagenesis
studies revealed the role of these amino acids in anion selectivity: Gu et al.
32 measured a weak cation selectivity when E111 and K147 were mutated into
neutral asparagines; Maglia et al. 44 found that the mutation of E111 into
asparagines yields a slight increase into anion selectivity; Mohammad and
Movileanu 33 observed a cation selectivity when K147 were mutated into
negatively charged aspartic acids.

We constructed three systems with modified charges: with protonated glutamic
acids E111 (E111N, resulting in a total top charge of $+7e$) with deprotonated
lysines K147 (K147N, top charge $-7e$), and with both amino acids neutralized
(E111N-K147N, null top charge as the WT-pore).

(A)0cm-3.9cm (B)0cm-3.9cm

Figure 6: $IV$ curves for $\alpha$HL with a modified constriction: E111N (A)
and K147N (B).

Figure 6 shows the $IV$ curves obtained for the E111N (A) and K147N (B)
systems. As observed on Fig. 6A and Table 1, E111N $\alpha$HL presents a
current asymmetry with a ratio $r=1.43$ and an anion selectivity similar to
the native protein ($s=0.71$), therefore suggesting that E111 has almost no
effect on selectivity or rectification. Yet, it should be noticed that, for
this system, the anchoring of the pore inside the membrane is weakened,
resulting in system crashes for high electric fields ($E_{z}=\pm 0.04$
V.nm-1). This result points out the important role of E111 in the stability of
the pore insertion inside the membrane. The total current for E111N is smaller
than the WT-protein, as shown by conductances in Table 1. This result is not
consistent with previous experiments of E111 mutations into asparagines for
which the total current slightly increased 44. However, this discrepancy can
be explained, as suggested by Stoddard et al. 7, by the fact that asparagines
are smaller residues than glutamic acids and the mutated pore diameter is
therefore a little wider than the WT at the constriction. The difference in
polarity between these residues can also explain this difference. Asparagines,
due to their polarity, do probably not allow the accumulation of anions at the
constriction as for the neutral glutamic acids (as explained below and visible
on Fig. 7A), resulting in a larger diameter and a higher current.

Figure 6B shows that the K147N pore presents a slightly weaker current
rectification than the WT-pore ($r=1.21$) and no ionic selectivity ($s=0.96)$.
Cationic and anionic currents are both asymmetric in this case
($r_{\text{Cl}}=1.29$), contrary to the WT pore for which the Cl- current is
nearly symmetric. Therefore, there is a rectification effect even without
selectivity. This result suggests that K147 plays a determinant role in the
ionic selectivity of $\alpha$HL. This is in accordance with a previous
mutagenesis study 33 which mutated K147 into aspartic acids and measured a
cation selectivity.

Finally, the E111N-K147N pore looses the ionic selectivity but also the
current asymmetry ($s=1.01$, $r=0.99$). The absence of rectification in this
case is due to both symmetric cation and anionic currents, whereas for the WT-
pore, rectification is due to the asymmetry in the Na+ current. It is in
accordance with a previous experimental study where E111 and K147 were mutated
to asparagines and current rectification disappeared 45. We can again notice
that anion selectivity and rectification effects are intertwined, and the
modification of one feature can influence the other.

The total current is greater for the two systems in which K147 was
deprotonated compared to the WT-pore, showing conductances close to the
neutralized pore (see Table 1). Since these systems present no anion
selectivity, the cationic current is greater than for the WT pore, yielding a
higher total current. This is consistent with several previous theoretical 22
and experimental studies 45, 46, but other experiments 32, 6, 7 did not
measure any increase of the ionic flow.

(A)0.2cm-2.0cm (B)0.2cm-2.5cm

Figure 7: Ionic density maps for the modified $\alpha$-hemolysins E111N (A)
and K147N (B) in the presence of a negative electric field $E_{z}=-0.03$ V/nm.

The average ionic density maps of the E111N (A) and K147N (B) systems, with a
negative electric field, are plotted on Figure 7. Both systems present a
current rectification like the WT-pore, but K147N show no anion selectivity,
whereas E111N is anion selective. Indeed, we observe a completely different
picture at the constriction between the two pores. The E111N (Fig. 7A) shows a
very high anion density around the constriction. This anion distribution
compensates the global positive charge of this region and prevents cations
from passing, hence responsible for the anion selectivity. On the contrary,
for the K147N, as shown on Fig. 7B, there are no anions at the constriction
site, and a high cation density is observed. The cations are allowed to pass
in this case, resulting in the absence of anion selectivity. In addition, the
channel is wider for K147N than for E111N, explaining that the total current
is greater than WT and E111N.

At the trans part, with a negative potential, the two maps reveal a positive
peak similar to the WT-pore. There are more cations located at the trans
extremity, which can explain that the rectification effect is still observed
for K147N, although we could expect that the presence of a cationic current
with negative electric fields would induce the absence of rectification (the
anionic current is much smaller for negative electric fields, as shown by the
$r_{\text{Cl}}=1.29$ ratio). Consequently, we remarked that charge
modifications at one stem extremity can influence the ionic density not only
close to the charges but also at the opposite extremity. For all non selective
pores, no anions were found at the constriction. Furthermore, with negative
electric fields, there is an increase in the number of cations at the trans
side (see Supplemental Information for the total number of ions in the stem of
K147N).

Our results confirm the role of the K147 charges in the anion selectivity of
$\alpha$HL. We also found that, like with modified charges at the trans part,
and according to the ionic density maps, it is the detailed residue charges of
the constriction part which are responsible for the anion selectivity, and not
the global charge of the restriction region.

Table 1: Conductance values from CG-MD simulations calculated from linear
regression on $IV$ curves. The total charge of the pore, the asymmetry ratio
$r=\frac{G^{+}}{G^{-}}$, the selectivity ratio
$s=\frac{G_{\text{Na}}^{+}}{G_{\text{Cl}}^{+}}$ and the asymmetry ratio for
the anionic current $r_{\text{Cl}}$ are indicated for each protein. The
selectivity ratio $s$ was evaluated using conductances for simulations with a
positive voltage (+).

| $Q_{\text{tot}}$ | $G^{+}$ | $G^{-}$ | $r$ | $s$ | $r_{\text{Cl}}$
---|---|---|---|---|---|---
| (e) | (pS) | (pS) | $G^{+}/G^{-}$ | $G_{\text{Na}}^{+}/G_{\text{Cl}}^{+}$ | $G_{\text{Cl}}^{+}/G_{\text{Cl}}^{-}$
WT | +7 | 68.3 | 50.7 | 1.35 | 0.60 | 0.91
Neutral* | 0 | 117.1 | 113.5 | 1.03 | 0.85 | 1.03
D127N* | +14 | 20.3 | 33.8 | 0.60 | 0.40 | 0.43
D128N* | +14 | 33.2 | 33.1 | 1.00 | 0.48 | 0.66
D127N-D128N-K131N* | +14 | 36.2 | 34.7 | 1.04 | 0.29 | 0.88
K131N* | 0 | 33.8 | 35.2 | 0.96 | 0.49 | 0.65
D127N-K131N* | +7 | 40.0 | 66.5 | 0.60 | 0.38 | 0.44
D128N-K131N* | +7 | 33.9 | 55.2 | 0.61 | 0.32 | 0.48
D127N-D128N* | +21 | 38.0 | 35.1 | 1.08 | 0.40 | 0.80
E111N-K147N | +7 | 120.0 | 119.4 | 1.00 | 0.99 | 1.01
E111N | +14 | 54.6 | 38.2 | 1.43 | 0.71 | 0.90
K147N | 0 | 110.3 | 91.5 | 1.21 | 0.96 | 1.29

* Simulations run with Position Restraints on the backbone beads

##  3 Conclusion

In this study, we characterized the role of certain amino acids of the
$\alpha$-hemolysin pore using coarse-grained molecular dynamics. By computing
$IV$-curves and ionic density maps of a very exhaustive panel of $\alpha$HL
charge modifications, we could reproduce experimentally observed behaviors as
well as all-atom MD simulations results. We demonstrated that the current
rectification is due to the trans side charges, in particular D128. On the
other hand, we showed that the anion selectivity arises from the K147 residues
located at the cis constriction. Even though these two effects are closely
linked, they are essentially due to subtle changes of the electrostatic
potential and not to the global charge of the pore or of the channel
extremities. The use of ionic density maps gives some microscopic insight of
these phenomenon confirming the importance of the local potential due to D128
and K147 in, respectively, current rectification and anion selectivity.

This study demonstrates the ability of CG-MD, with a polarizable solvent, to
reliably model nanopore currents using voltages closer to experimental
conditions and smaller than AA simulations. It can describe the molecular
insight of processes with a timescale getting closer to translocation
experiments and it allows the study of numerous mutations. This modeling
system could be used to model nanopore molecular transport such as DNA
translocation or could be applied to other nanopores like aerolysin.

##  4 Methods

###  4.1 Simulations setup

Different systems were studied, either with a native pore or with a modified
$\alpha$-hemolysin. They were all prepared with the same protocol whose first
steps are identical to our previous paper 34.

First, the crystallographic structure of $\alpha$-hemolysin was taken from the
Protein Data Bank (PDB entry: 7AHL) 47 The coordinates of the missing atoms
were recovered using the _pdb2pdr_ software 48 and the water molecules were
removed. Using the _propka_ software 49, we determined the protonation states
of charged residues: all histidines were therefore considered neutral.

The protein was reduced in coarse grains using the _martinize.py_ script
provided by MARTINI, with the Martini 2.2 force field 50. The ElNeDyn approach
from MARTINI was applied on the protein 35, adding an elastic network on the
backbone and between the different chains of the $\alpha$-hemolysin, which
consists of springs with 500 kJ.mol-1.nm-2 elastic bond strength connecting
grains within a 0.9 nm cutoff from each other.

After a short minimization, the backbone of the resulting protein structure
was fitted on the one of a coarse-grained $\alpha$-hemolysin inserted in a
dipalmitoylphosphatidylcholine (DPPC) bilayer obtained by Pr. Sansom 51. This
last protein structure was then deleted resulting in a system composed of our
$\alpha$-hemolysin structure inserted in a DPPC bilayer of 756 molecules in a
simulation box initially set at 15.85 x 15.85 x 20 nm. A short minimization
was performed and the system was solvated with 28,758 W MARTINI CG water
molecules using the _genbox_ option from the Gromacs package (version 4.6).
Another minimization was then carried out followed by two equilibration steps:
one of 1 ns in the NVT ensemble (using a 10 fs timestep) and the other of 5 ns
in the NPT ensemble (using a 20 fs timestep).

Afterwards, the W MARTINI water was transformed into polarizable water PW 52
using the _triple_w.py_ script provided by MARTINI, and the system underwent a
short minimization. The charges of the amino acids of the protein were then
modified as required and we performed two equilibration steps (using Gromacs
package version 5.1.4): 2 ns at NVT and 10 ns at NPT, both with position
restraints on the protein backbone. Next, we added the ions with the _genion_
tool of Gromacs in order to obtain a neutral system with 1 M ionic
concentration. For example, for a neutral $\alpha$-hemolysin, we added 2054 of
both Na+ and Cl- MARTINI ions. We carried out another equilibration step of 10
ns at NPT with position restraints on the protein backbone. The resulting
systems were used as initial state for the simulations.

It should be noted that the MARTINI ions named Na+ and Cl- are totally
symmetrical in terms of mass ($m=72$ g.mol-1) and size, but carry opposite
charges. We can therefore compare our coarse-grained simulations to
experiments or all-atom simulations with KCl ions.

###  4.2 MD methods

All our molecular dynamics simulations were performed using the GROMACS
package programs (version 5.1.4) with the MARTINI force field
(martini-v.2.2P). Periodic boundaries conditions were applied on the system.
We used Particle-mesh Ewald (PME) method for long-range electrostatics
calculations, with a 2 Å Fourier grid spacing and a direct space cut-off
radius of 13 Å. Simulations were performed on the NPT ensemble, at a
temperature of 320 K and 1 bar pressure. Nosé-Hoover thermostat and Parinello-
Rahman barostat were applied. A timestep of 20 fs was used for all the
simulations.

Different electric fields $E_{z}$ were applied in the $z$ direction,
perpendicularly to the membrane plan, for each system, going from $-0.04$
V.nm-1 to $+0.04$ V.nm-1: 0, $\pm 0.01$, $\pm 0.02$, $\pm 0.03$ and $\pm 0.04$
V.nm-1. Each simulation was conducted for 1.5 $\mu$s with a constant electric
field applied.

###  4.3 Ionic current analyses

For each applied electric field, we counted the ions crossing a horizontal
plan at the trans side of the pore which allowed us to obtain cumulative
current over time, similarly to Refs 53, 34. The average current intensity was
calculated from the slopes of the linear regressions on the cumulative current
for 250 ns intervals with a shift of 50 ns. We computed this current intensity
on 500 ns from 1 to 1.5 $\mu$s of simulation. The applied voltage
corresponding to each applied electric field was evaluated using $\Delta
V={E_{z}}\times{L_{z}}$, with $\Delta V$ the voltage in V, $E_{z}$ the applied
electric field in V.nm-1 and $L_{z}$ the average length of the box in $z$
direction, during the relevant time window, in nm. From these calculations of
average current intensities and applied voltages, we plotted the
characteristic $IV$ curves of the pores, corresponding to the current
intensities as a function of the voltage applied on the system.

Intensities generated by either the Na+ or the Cl- flows were also obtained as
a function of the applied electric fields, by the same method, and were
represented separately on the $IV$ curves.

The positive ($G^{+}$) and negative ($G^{-}$) conductances were computed using
the slope of linear regressions on these $IV$ curves for the positive or
negative voltages, respectively.

The current asymmetry was calculated using the ratio $r=\frac{G^{+}}{G^{-}}$
and the ionic selectivity was deduced from the ratio
$s=\frac{{G_{Na^{+}}}^{+}}{{G_{Cl^{-}}}^{+}}$ with ${G_{Na^{+}}}^{+}$ and
${G_{Cl^{-}}}^{+}$ the positive conductances for each ion species.

###  4.4 Ionic charge densities analyses

We used _densmap_ tool from Gromacs package to plot the average ionic
densities of the different systems from 1 to 1.5 µs of simulation. This tool
can compute the average particle density during the simulation around an axis
resulting in 2D-maps. Therefore, we defined the axis by using the residues at
the bottom (D127, D128, T129, G130 and K131) and at the top (N17) of the pore.
Moreover, we determined Na+ and Cl- densities with a bin size of 0.25 nm and
with the mirror option of _densmap_ to represent the density on both sides of
the axis, leading to in a symmetric density centered on the pore axis.

We calculated the charge density by combining Na+ and Cl- individual densities
such as: $d=q_{Na}\times d_{Na}+q_{Cl}\times d_{Cl}$ with $q$ the ion charge
and $d$ the density in one bin. Therefore, when there are more cations than
anions in one bin, the charge density is positive and conversely for negative
density.

This work was granted access to the HPC resources of CINES under the
allocations A0020707139 (2017), A0040707139 (2018) and A0060707139 (2019) made
by GENCI.

Number of ions in function of applied voltage on WT and K147N $\alpha$HL
pores. $IV$ curves for D127ND128N-K131N*, D127N*, K131N*, D127N-D128N*,
D128N-K131N* and E111N-K147N $\alpha$HL. Ionic density maps for the neutral
constrained pore and for D127N-D128N-K131N*, E111N-K147N, K131N*, D127N-D128N*
and D128N-K131N* $\alpha$HL.

## References

  * Kasianowicz et al. 1996 Kasianowicz, J. J.; Brandin, E.; Branton, D.; Deamer, D. W. Characterization of Individual Polynucleotide Molecules Using a Membrane Channel. _Proc. Natl. Acad. Sci. U. S. A._ 1996, _93_ , 13770–13773
  * Dudko et al. 2010 Dudko, O. K.; Mathé, J.; Meller, A. Nanopore Force Spectroscopy Tools for Analyzing Single Biomolecular Complexes. _Methods Enzymol._ 2010, _475_ , 565–589
  * Meller et al. 2000 Meller, A.; Nivon, L.; Brandin, E.; Golovchenko, J.; Branton, D. Rapid Nanopore Discrimination Between Single Polynucleotide Molecules. _Proc. Natl. Acad. Sci. U. S. A._ 2000, _97_ , 1079–1084
  * Oukhaled et al. 2007 Oukhaled, G.; Mathé, J.; Biance, A.-L.; Bacri, L.; Betton, J.-M.; Lairez, D.; Pelta, J.; Auvray, L. Unfolding of Proteins and Long Transient Conformations Detected by Single Nanopore Recording. _Phys. Rev. Lett._ 2007, _98_
  * Oukhaled et al. 2008 Oukhaled, G.; Bacri, L.; Mathé, J.; Pelta, J.; Auvray, L. Effect of Screening on the Transport of Polyelectrolytes through Nanopores. _EPL_ 2008, _82_ , 48003
  * Stoddart et al. 2009 Stoddart, D.; Heron, A. J.; Mikhailova, E.; Maglia, G.; Bayley, H. Single-Nucleotide Discrimination in Immobilized DNA Oligonucleotides with a Biological Nanopore. _Proc. Natl. Acad. Sci. U. S. A._ 2009, _106_ , 7702–7707
  * Stoddart et al. 2010 Stoddart, D.; Heron, A. J.; Klingelhoefer, J.; Mikhailova, E.; Maglia, G.; Bayley, H. Nucleobase Recognition in ssDNA at the Central Constriction of the $\alpha$-Hemolysin Pore. _Nano Lett._ 2010, _10_ , 3633–3637
  * Deamer et al. 2016 Deamer, D.; Akeson, M.; Branton, D. Three Decades of Nanopore Sequencing. _Nat. Biotechnol._ 2016, _34_ , 518
  * Mathé et al. 2004 Mathé, J.; Visram, H.; Viasnoff, V.; Rabin, Y.; Meller, A. Nanopore Unzipping of Individual DNA Hairpin Molecules. _Biophys. J._ 2004, _87_ , 3205–3212
  * Mathé et al. 2005 Mathé, J.; Aksimentiev, A.; Nelson, D. R.; Schulten, K.; Meller, A. Orientation Discrimination of Single-Stranded DNA Inside the $\alpha$-Hemolysin Membrane Channel. _Proc. Natl. Acad. Sci. U. S. A._ 2005, _102_ , 12377–12382
  * Muzard et al. 2010 Muzard, J.; Martinho, M.; Mathé, J.; Bockelmann, U.; Viasnoff, V. DNA Translocation and Unzipping through a Nanopore: Some Geometrical Effects. _Biophys. J._ 2010, _98_ , 2170–2178
  * Di Marino et al. 2015 Di Marino, D.; Bonome, E. L.; Tramontano, A.; Chinappi, M. All-Atom Molecular Dynamics Simulation of Protein Translocation through an $\alpha$-Hemolysin Nanopore. _J. Phys. Chem. Lett._ 2015, _6_ , 2963–2968
  * Di Muccio et al. 2019 Di Muccio, G.; Rossini, A. E.; Di Marino, D.; Zollo, G.; Chinappi, M. Insights into Protein Sequencing with an $\alpha$-Hemolysin Nanopore by Atomistic Simulations. _Sci. Rep._ 2019, _9_
  * Bonome et al. 2019 Bonome, E. L.; Cecconi, F.; Chinappi, M. Translocation Intermediates of Ubiquitin through an $\alpha$-Hemolysin Nanopore: Implications for Detection of Post-Translational Modifications. _Nanoscale_ 2019, _11_ , 9920–9930
  * Fennouri et al. 2018 Fennouri, A.; Ramiandrisoa, J.; Bacri, L.; Mathé, R., Jérôme ans Danie Comparative Biosensing of Glycosaminoglycan Hyaluronic Acid Oligo- and Polysaccharides Using Aerolysin and $\alpha$-Hemolysin Nanopores. _Eur. Phys. J. E_ 2018, _41_ , 1–7
  * Krasilnikov et al. 2006 Krasilnikov, O. V.; Rodrigues, C. G.; Bezrukov, S. M. Single Polymer Molecules in a Protein Nanopore in the Limit of a Strong Polymer-Pore Attraction. _Phys. Rev. Lett._ 2006, _97_ , 018301
  * Reiner et al. 2010 Reiner, J. E.; Kasianowicz, J. J.; Nablo, B. J.; Robertson, J. W. F. Theory for Polymer Analysis Using Nanopore-Based Single-Molecule Mass Spectrometry. _Proc. Natl. Acad. Sci. U. S. A._ 2010, _107_ , 12080–12085
  * Merzlyak et al. 2005 Merzlyak, P. G.; Capistrano, M.-F. P.; Valeva, A.; Kasianowicz, J. J.; Krasilnikov, O. V. Conductance and Ion Selectivity of a Mesoscopic Protein Nanopore Probed with Cysteine Scanning Mutagenesis. _Biophys. J._ 2005, _89_ , 3059–3070
  * Piguet et al. 2014 Piguet, F.; Discala, F.; Breton, M.-F.; Pelta, J.; Bacri, L.; Oukhaled, A. Electroosmosis through $\alpha$-Hemolysin That Depends on Alkali Cation Type. _J. Phys. Chem. Lett._ 2014, _5_ , 4362–4367
  * Payet et al. 2015 Payet, L.; Martinho, M.; Merstorf, C.; Pastoriza-Gallego, M.; Pelta, J.; Viasnoff, V.; Auvray, L.; Muthukumar, M.; Mathé, J. Temperature Effect on Ionic Current and ssDNA Transport through Nanopores. _Biophys. J._ 2015, _109_ , 1600–1607
  * Misakian and Kasianowicz 2003 Misakian, M.; Kasianowicz, J. J. Electrostatic Influence on Ion Transport through the aHL Channel. _J. Membr. Biol._ 2003, _195_ , 137–146
  * Noskov et al. 2004 Noskov, S. Y.; Im, W.; Roux, B. Ion Permeation through the $\alpha$-Hemolysin Channel: Theoretical Studies Based on Brownian Dynamics and Poisson-Nernst-Plank Electrodiffusion Theory. _Biophys. J._ 2004, _87_ , 2299–2309
  * Aksimentiev and Schulten 2005 Aksimentiev, A.; Schulten, K. Imaging $\alpha$-Hemolysin with Molecular Dynamics: Ionic Conductance, Osmotic Permeability, and the Electrostatic Potential Map. _Biophys. J._ 2005, _88_ , 3745–3761
  * Millar et al. 2008 Millar, C.; Madathil, R.; Beckstein, O.; Sansom, M. S. P.; Roy, S.; Asenov, A. Brownian Simulation of Charge Transport in $\alpha$-Haemolysin. _J. Comput. Electron._ 2008, _7_ , 28–33
  * Bhattacharya et al. 2011 Bhattacharya, S.; Muzard, J.; Payet, L.; Mathé, J.; Bockelmann, U.; Aksimentiev, A.; Viasnoff, V. Rectification of the Current in $\alpha$-Hemolysin Pore Depends on the Cation Type: The Alkali Series Probed by Molecular Dynamics Simulations and Experiments. _J. Phys. Chem. C_ 2011, _115_ , 4255–4264
  * Gamble et al. 2014 Gamble, T.; Decker, K.; Plett, T. S.; Pevarnik, M.; Pietschmann, J.-F.; Vlassiouk, I.; Aksimentiev, A.; Siwy, Z. S. Rectification of Ion Current in Nanopores Depends on the Type of Monovalent Cations: Experiments and Modeling. _J. Phys. Chem. C_ 2014, _118_ , 9809–9819
  * Manara et al. 2015 Manara, R. M. A.; Guy, A. T.; Wallace, E. J.; Khalid, S. Free-Energy Calculations Reveal the Subtle Differences in the Interactions of DNA Bases with $\alpha$-Hemolysin. _J. Chem. Theory Comput._ 2015, _11_ , 810–816, PMID: 26579606
  * Bonome et al. 2017 Bonome, E. L.; Cecconi, F.; Chinappi, M. Electroosmotic Flow through an $\alpha$-Hemolysin Nanopore. _Microfluid. Nanofluid._ 2017, _21_
  * Zhou et al. 2020 Zhou, W.; Qiu, H.; Guo, Y.; Guo, W. Molecular Insights into Distinct Detection Properties of $\alpha$-Hemolysin, MspA, CsgG, and Aerolysin Nanopore Sensors. _J. Phys. Chem. B_ 2020, _124_ , 1611–1618, PMID: 32027510
  * Menestrina 1986 Menestrina, G. Ionic Channels Formed by Staphylococcus aureus Alpha-Toxin: Voltage-Dependent Inhibition by Divalent and Trivalent Cations. _J. Membr. Biol._ 1986, _90_ , 177–190
  * Bonthuis et al. 2006 Bonthuis, D. J.; Zhang, J.; Hornblower, B.; Mathé, J.; Shklovskii, B. I.; Meller, A. Self-Energy-Limited Ion Transport in Subnanometer Channels. _Phys. Rev. Lett._ 2006, _97_ , 128104
  * Gu et al. 2000 Gu, L.-Q.; Dalla Serra, M.; Vincent, J. B.; Vigh, G.; Cheley, S.; Braha, O.; Bayley, H. Reversal of Charge Selectivity in Transmembrane Protein Pores by Using Noncovalent Molecular Adapters. _Proc. Natl. Acad. Sci. U. S. A._ 2000, _97_ , 3959–3964
  * Mohammad and Movileanu 2010 Mohammad, M. M.; Movileanu, L. Impact of Distant Charge Reversals within a Robust $\beta$-Barrel Protein Pore. _J. Phys. Chem. B_ 2010, _114_ , 8750–8759
  * Basdevant et al. 2019 Basdevant, N.; Dessaux, D.; Ramirez, R. Ionic Transport through a Protein Nanopore: a Coarse-Grained Molecular Dynamics Study. _Sci. Rep._ 2019, _9_ , 15740
  * Periole et al. 2009 Periole, X.; Cavalli, M.; Marrink, S.-J.; Ceruso, M. A. Combining an Elastic Network With a Coarse-Grained Molecular Force Field: Structure, Dynamics, and Intermolecular Recognition. _J. Chem. Theory Comput._ 2009, _5_ , 2531–2543
  * Gumbart et al. 2012 Gumbart, J.; Khalili-Araghi, F.; Sotomayor, M.; Roux, B. Constant Electric Field Simulations of the Membrane Potential Illustrated with Simple Systems. _Biochim. Biophys. Acta, Biomembr._ 2012, _1818_ , 294–302
  * Modi et al. 2012 Modi, N.; Winterhalter, M.; Kleinekathöfer, U. Computational Modeling of Ion Transport through Nanopores. _Nanoscale_ 2012, _4_ , 6166–6180
  * Aksimentiev 2010 Aksimentiev, A. Deciphering Ionic Current Signatures of DNA Transport through a Nanopore. _Nanoscale_ 2010, _2_ , 468–483
  * Crozier et al. 2001 Crozier, P. S.; Henderson, D.; Rowley, R. L.; Busath, D. D. Model Channel Ion Currents in NaCl-Extended Simple Point Charge Water Solution with Applied-Field Molecular Dynamics. _Biophys. J._ 2001, _81_ , 3077–3089
  * Walker et al. 1992 Walker, B.; Krishnasastry, M.; Zorn, L.; Kasianowicz, J.; Bayley, H. Functional Expression of the $\alpha$-Hemolysin of Staphylococcus aureus in Intact Escherichia coli and in Cell Lysates. Deletion of Five C-Terminal Amino Acids Selectively Impairs Hemolytic Activity. _J. Biol. Chem._ 1992, _267_ , 10902–10909
  * Simakov and Kurnikova 2010 Simakov, N. A.; Kurnikova, M. G. Soft Wall Ion Channel in Continuum Representation with Application to Modeling Ion Currents in $\alpha$-Hemolysin. _J. Phys. Chem. B_ 2010, _114_ , 15180–15190
  * Dyrka et al. 2013 Dyrka, W.; Bartuzel, M. M.; Kotulska, M. Optimization of 3D Poisson-Nernst-Planck Model for Fast Evaluation of Diverse Protein Channels: 3D PNP for Diverse Protein Channels. _Proteins: Struct., Funct., Bioinf._ 2013, _81_ , 1802–1822
  * Boukhet et al. 2016 Boukhet, M.; Piguet, F.; Ouldali, H.; Pastoriza-Gallego, M.; Pelta, J.; Oukhaled, A. Probing Driving Forces in Aerolysin and $\alpha$-Hemolysin Biological Nanopores: Electrophoresis Versus Electroosmosis. _Nanoscale_ 2016, _8_ , 18352–18359
  * Maglia et al. 2008 Maglia, G.; Restrepo, M. R.; Mikhailova, E.; Bayley, H. Enhanced Translocation of Single DNA Molecules through $\alpha$-Hemolysin Nanopores by Manipulation of Internal Charge. _Proc. Natl. Acad. Sci. U. S. A._ 2008, _105_ , 19720–19725
  * Rincon-Restrepo et al. 2011 Rincon-Restrepo, M.; Mikhailova, E.; Bayley, H.; Maglia, G. Controlled Translocation of Individual DNA Molecules through Protein Nanopores with Engineered Molecular Brakes. _Nano Lett._ 2011, _11_ , 746–750
  * Ayub and Bayley 2012 Ayub, M.; Bayley, H. Individual RNA Base Recognition in Immobilized Oligonucleotides Using a Protein Nanopore. _Nano Lett._ 2012, _12_ , 5637–5643
  * Berman et al. 2000 Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. _Nucleic Acids Res._ 2000, _28_ , 235–242
  * Dolinsky et al. 2004 Dolinsky, T. J.; Nielsen, J. E.; McCammon, J. A.; Baker, N. A. PDB2PQR: an Automated Pipeline for the Setup of Poisson-Boltzmann Electrostatics Calculations. _Nucleic Acids Res._ 2004, _32_ , W665–W667
  * Olsson et al. 2011 Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. _J. Chem. Theory Comput._ 2011, _7_ , 525–537
  * de Jong et al. 2013 de Jong, D. H.; Singh, G.; Bennett, W. F. D.; Arnarez, C.; Wassenaar, T. A.; Schäfer, L. V.; Periole, X.; Tieleman, D. P.; Marrink, S. J. Improved Parameters for the Martini Coarse-Grained Protein Force Field. _J. Chem. Theory Comput._ 2013, _9_ , 687–697
  * Scott et al. 2008 Scott, K. A.; Bond, P. J.; Ivetac, A.; Chetwynd, A. P.; Khalid, S.; Sansom, M. S. P. Coarse-Grained MD Simulations of Membrane Protein-Bilayer Self-Assembly. _Structure_ 2008, _16_ , 621–630
  * Yesylevskyy et al. 2010 Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable Water Model for the Coarse-Grained MARTINI Force Field. _PLoS Comput. Biol._ 2010, _6_ , e1000810
  * Guros et al. 2018 Guros, N. B.; Balijepalli, A.; Klauda, J. B. The Role of Lipid Interactions in Simulations of the $\alpha$-Hemolysin Ion-Channel-Forming Toxin. _Biophys. J._ 2018, _115_ , 1720–1730


# Time-dependent study of $K_{S}\to\pi^{+}\pi^{-}$ decays for flavour physics
measurements

P. Pakhlov<EMAIL_ADDRESS>P.N. Lebedev Physical Institute of the RAS,
Moscow, Russia  Higher School of Economics (National Research University),
Moscow, Russia  V. Popov<EMAIL_ADDRESS>P.N. Lebedev Physical Institute
of the RAS, Moscow, Russia

###### Abstract

Nowadays High Energy Physics experiments can accumulate unprecedented
statistics of heavy flavour decays that allows to apply new methods, based on
the study of very rare phenomena, which used to be just desperate. In this
paper we propose a new method to measure composition of
$K^{0}$-$\overline{K}{}^{0}$, produced in a decay of heavy hadrons. This
composition contains important information, in particular about weak and
strong phases between amplitudes of the produced $K^{0}$ and
$\overline{K}{}^{0}$. We consider possibility to measure these parameters with
time-dependent $K^{0}\to\pi^{+}\pi^{-}$ analysis. Due to $CP$-violation in
kaon mixing time-dependent decay rates of $K^{0}$ and $\overline{K}{}^{0}$
differ, and the initial amplitudes revealed in the $CP$-violating decay
pattern. We perform phenomenological study of $K^{0}$ decay evolution
initially produced as a combination
$a\left|K^{0}(t)\right\rangle+b\left|\overline{K}{}^{0}(t)\right\rangle$,
where $a$ and $b$, complex amplitudes, could also be dependent on decay time
of heavy mother particle. In particular we consider cases of charmed hadrons
decays: $D^{+}\to K^{0}\pi^{+}$, $D_{s}^{+}\to K^{0}K^{+}$, $\Lambda\to
pK^{0}$ and with some assumptions $D^{0}\to K^{0}\pi^{0}$. This can be used to
test the sum rule for charmed mesons and to obtain input for the full
constraint of the two body amplitudes of $D$-mesons.

##  I Introduction

The final states of heavy hadron decays sometimes differ only by the
replacement of the strange $K^{0}$ meson by its antiparticle due to the
contribution of different diagrams or neutral meson mixing, thus representing
a superposition of the strange and antistrange states. An example is the decay
$D^{+}\to\overline{K}{}^{0}(K^{0})\pi^{+}$, where both $\overline{K}{}^{0}$
and $K^{0}$ are produced due to the presence of both Cabibbo favourite (CF)
and doubly Cabibbo suppressed (DCS) amplitudes.

Considering charm hadron decays from the point of view of flavour $SU(3)_{f}$
symmetry one can obtain the following sum rules Grossman:2012ry :

$\displaystyle\sqrt{2}A_{D^{0}\to\overline{K}{}^{0}\pi^{0}}+A_{D^{0}\to
K^{-}\pi^{+}}-A_{D^{+}\to\overline{K}{}^{0}\pi^{+}}=0,$ (1)
$\displaystyle\sqrt{2}A_{D^{0}\to K^{0}\pi^{0}}+A_{D^{0}\to
K^{+}\pi^{-}}+\sqrt{2}A_{D^{+}\to K^{+}\pi^{0}}-A_{D^{+}\to K^{0}\pi^{+}}=0.$
(2)

These relations are the isospin sum rules and both are broken at the same
level of $\mathcal{O}\left((m_{u}-m_{d})/\Lambda_{QCD}\right)\sim 1\%$.
Amplitudes involved in eq. (1, 2) are of the same order of Cabibbo suppression
and correspond to CF and DCS decay amplitudes, respectively. These sum rules
could be illustrated as shown in Fig. 1.

$\mathbf{\sqrt{2}A_{D^{0}\to\overline{K}^{0}\pi^{0}}}$$\mathbf{A_{D^{0}\to
K^{-}\pi^{+}}}$$\mathbf{A_{D^{+}\to\overline{K}^{0}\pi^{+}}}$$\mathbf{\sqrt{2}A_{D^{0}\to
K^{0}\pi^{0}}}$$\mathbf{A_{D^{0}\to
K^{+}\pi^{-}}}$$\mathbf{\sqrt{2}A_{D^{+}\to
K^{+}\pi^{0}}}$$\mathbf{A_{D^{+}\to
K^{0}\pi^{+}}}$$\delta^{00}$$\delta^{+-}$$\delta^{0+}$$\mathrm{(a)}$$\mathrm{(b)}$
Figure 1: $SU_{f}(3)$ sum rules for CF (a) and DCS (b) decay amplitudes.

There is another sum based on $SU(3)_{f}$ that is particularly interesting
since here CF and DCS amplitudes are mixed together:

$A_{D^{+}\to\overline{K}{}^{0}\pi^{+}}-A_{D_{s}^{+}\to\overline{K}{}^{0}K^{+}}-\frac{A_{D^{+}\to\overline{K}{}^{0}K^{+}}}{\lambda}+\frac{A_{D_{s}^{+}\to\overline{K}{}^{0}\pi^{+}}}{\lambda}+\frac{A_{D^{+}\to
K^{0}\pi^{+}}}{\lambda^{2}}-\frac{A_{D_{s}^{+}\to
K^{0}K^{+}}}{\lambda^{2}}=0,$ (3)

where $\lambda=\sin{\theta_{12}}$ – is the CKM parameter and $\theta_{12}$ is
the Cabibbo angle.

Recently, the question of validity of $SU(3)_{f}$ rules and accuracy of their
approximation has become urgent as they are widely used to explain the
anomalously large $CP$ violation observed in $D^{0}$ decays by LHCb experiment
Aaij:2019kcg :

$\Delta A_{CP}=A_{CP}(K^{+}K^{-})-A_{CP}(\pi^{+}\pi^{-})=(-15.4\pm 2.9)\times
10^{-4},$ (4)

Many papers trying to reconcile LHCb result with the Standard Model suppose
that such discrepancy could be explained by the enhanced penguin amplitudes
Golden:1989qx ; Hiller:2012xm ; Brod:2012ud ; Franco:2012ck . In these studies
it was argued that $SU(3)_{f}$ breaking not necessary happens at the amplitude
level but rather could be explained by FSI generated phases.

Structure of charm decay amplitudes could be described in terms of Wick
contractions of operators of effective Hamiltonian. In this work we are
particularly interested in amplitudes involving $K^{0}$ in the final state. In
the following we adopt results for amplitude expressions obtained in
Buccella:2019kpn , for Cabibbo-favourite amplitudes:

$\displaystyle A_{D^{0}\to K^{-}\pi^{+}}$ $\displaystyle=$
$\displaystyle\frac{1}{5}(3T-2C-K)e^{i\delta_{1/2}}+\frac{2}{5}(T+C+\kappa)$
$\displaystyle A_{D^{0}\to\overline{K}{}^{0}\pi^{0}}$ $\displaystyle=$
$\displaystyle-\frac{1}{5\sqrt{2}}(3T-2C-K)e^{i\delta_{1/2}}+\frac{3}{5\sqrt{2}}(T+C+\kappa)$
$\displaystyle A_{D^{+}\to\overline{K}{}^{0}\pi^{0}}$ $\displaystyle=$
$\displaystyle(T+C+\kappa)$ $\displaystyle
A_{D_{s}^{+}\to\overline{K}{}^{0}K^{+}}$ $\displaystyle=$
$\displaystyle-\frac{1}{5}(2T-3C+\Delta)e^{i\delta^{\prime}_{1}}+\frac{2}{5}(T+C+\kappa),$
(5)

and for doubly Cabibbo-suppressed amplitudes:

$\displaystyle A_{D^{0}\to K^{+}\pi^{-}}$ $\displaystyle=$
$\displaystyle-\frac{1}{5}(3T-2C+K)e^{i\delta_{1/2}}-\frac{2}{5}(T+C+\kappa^{\prime})$
$\displaystyle A_{D^{0}\to K^{0}\pi^{0}}$ $\displaystyle=$
$\displaystyle\frac{1}{5\sqrt{2}}(3T-2C+K)e^{i\delta_{1/2}}-\frac{3}{5\sqrt{2}}(T+C+\kappa^{\prime})$
$\displaystyle A_{D^{+}\to K^{0}\pi^{+}}$ $\displaystyle=$
$\displaystyle\frac{1}{5}(2T-3C+\Delta-K^{\prime})e^{i\delta_{1/2}}-\frac{2}{5}(T+C+\kappa^{\prime})$
$\displaystyle A_{D_{s}^{+}\to K^{0}K^{+}}$ $\displaystyle=$
$\displaystyle-(T+C+\kappa^{\prime}).$ (6)

Here $T$ and $C$ correspond to “tree”-level color-connected and color-
suppressed amplitudes, $K$ and $K^{\prime}$ are parameters corresponding to
non-conservation of strangeness changing currents, $\kappa$ and
$\kappa^{\prime}$ – parameters allowing $SU(3)_{f}$-breaking in CF and DCS
amplitudes. Also two phases are present $\delta_{1/2}$, $\delta_{1}^{\prime}$
corresponding to $I=\frac{1}{2}$ and $I=1$ amplitudes respectively.

Buccella et al. Buccella:2019kpn  performed fit to the above mentioned
amplitude parameters based on observed values of $CP$ asymmetries and
branching ratios for charm hadron decays, and made predictions for the strong
phase difference in the $D^{0}\to K^{\pm}\pi^{\mp}$ decay channel:

$\delta^{+-}=(3.14\pm 5.69)^{\circ}.$ (7)

Furthermore, this study can be extended to obtain predictions for other strong
phase differences:

$\displaystyle\delta^{00}$ $\displaystyle=$ $\displaystyle(-3\pm 6)^{\circ};$
(8) $\displaystyle\delta^{0+}$ $\displaystyle=$ $\displaystyle(-76\pm
4)^{\circ};$ (9) $\displaystyle\delta_{s}^{0+}$ $\displaystyle=$
$\displaystyle(108\pm 4)^{\circ}.$ (10)

In this paper we propose a method that allows to measure strong phase
differences and amplitude ratios for final states with $K^{0}$-meson, hence
probe the validity of approaches used to explain LHCb result and in general
the sum rules. Such measurements will possibly allow to identify the source
and scale of $SU(3)_{f}$ breaking in charm hadron decays.

In our previous paper Pakhlov:2019wkn  we suggested to use semileptonic
$K^{0}$ decays and study their time evolution to measure the complex phase
between $K^{0}$ and $\overline{K}{}^{0}$ at the production point. The method
is based on the disentanglement of the production phase from the (known)
$K^{0}-\overline{K}{}^{0}$ mixing phase, using the later as a reference. While
the method can work properly, its experimental application is a challenge due
to unobserved neutrino in kaon decay. We realized that the similar sensitivity
could be achieved using the standard $K^{0}\to\pi^{+}\pi^{-}$ decays, which
are easy to reconstruct. In this paper we show that time evolution of the
state $a\left|K^{0}(t)\right\rangle+b\left|\overline{K}{}^{0}(t)\right\rangle$
decaying into $\pi^{+}\pi^{-}$ allows to extract the complex parameters $a$
and $b$ as well. In this case we utilize the $CP$-violating phase $\phi_{+-}$
as a reference for the proposed measurements. Both methods are based on the
rare neutral kaon effects: rare decays (semileptonic for short lived kaon
component) or rare interference (between short and long lived components).

##  II Method

Time evolution of the neutral kaon system could be described by Shrödinger
equation

$i{\partial}_{t}{K^{0}(t)\choose\overline{K}{}^{0}(t)}=\left({{\mathbf{M}}-\frac{i}{2}{\mathbf{\Gamma}}}\right){K^{0}(t)\choose\overline{K}{}^{0}(t)}\,,$
(11)

where effective Hamiltonian is a sum of absorptive and dispersive parts, $\bf
M$ and $\bf\Gamma$ are $2\times 2$ hermitian matrices. The Hamiltonian
eigenvalues could be written as follows:

$\displaystyle\lambda_{S,L}\equiv
m_{S,\,L}-i\frac{\Gamma_{S,\,L}}{2}=\left(M_{11}-i\frac{\Gamma_{11}}{2}\right)\pm\left(\frac{p}{q}\right)_{\\!\\!K}\left(M_{12}-i\frac{\Gamma_{12}}{2}\right)\,,$
(12)

where $m_{S,\,L}$, $\Gamma_{S,\,L}$ are masses and widths of the Hamiltonian
eigenstates $K^{0}_{S}$ and $K^{0}_{L}$, and parameters $p$, $q$ correspond to
the flavour admixtures of eigenstates defined by

$\left(\frac{p}{q}\right)_{\\!\\!K}^{2}=\frac{M_{12}-\frac{i}{2}\Gamma_{12}}{M_{12}^{*}-\frac{i}{2}\Gamma_{12}^{*}}\,.$
(13)

Since we consider only $\pi^{+}\pi^{-}$ final state it is convenient to use
$p/q=(1+\varepsilon)/(1-\varepsilon)$, where $\varepsilon$ describes the
$CP$-even component in $K^{0}_{L}$. Then amplitudes describing evolution of
initially pure flavour eigenstates in terms of $K^{0}_{S}$/$K^{0}_{L}$ could
be written as

$\displaystyle\left|K^{0}(t)\right\rangle$ $\displaystyle=$
$\displaystyle\frac{(1-\varepsilon)}{\sqrt{2}}\left[e^{-i\lambda_{S}t}\left|K_{S}\right\rangle+e^{-i\lambda_{L}t}\left|K_{L}\right\rangle\right],$
(14) $\displaystyle\left|\overline{K}{}^{0}(t)\right\rangle$ $\displaystyle=$
$\displaystyle\frac{(1+\varepsilon)}{\sqrt{2}}\left[e^{-i\lambda_{S}t}\left|K_{S}\right\rangle-e^{-i\lambda_{L}t}\left|K_{L}\right\rangle\right].$
(15)

Using these equations one can obtain for time-dependent decay rates

$\mathcal{R}(t)=\frac{1\mp
2\mathrm{Re}(\varepsilon)}{2}|A_{K_{S}\to\pi\pi}|^{2}\Big{[}e^{-\Gamma_{S}t}+|\eta_{+-}|^{2}e^{-\Gamma_{L}t}\\\
\pm 2|\eta_{+-}|e^{-\frac{1}{2}(\Gamma_{L}+\Gamma_{S})t}\cos{(\Delta
mt-\phi_{+-})}\Big{]},$ (16)

where the upper (lower) sign corresponds to initial pure $K^{0}$
($\overline{K}{}^{0}$), $\Delta m=m_{L}-m_{S}$ is a mass difference and we
introduced for the amplitude ratio parameter

$\frac{\left\langle\pi^{+}\pi^{-}\left|H\right|K^{0}_{L}\right\rangle}{\left\langle\pi^{+}\pi^{-}\left|H\right|K^{0}_{S}\right\rangle}=\eta_{+-}=|\eta_{+-}|e^{i\phi_{+-}}.$
(17)

The third interference term in eq. (16) basically allows to distinguish the
initial flavour of neutral kaon. Parameters $\eta_{+-}$ and $\phi_{+-}$ have
been measured with great precision and current world averages (assuming CPT
invariance) are Zyla:2020zbs : $\eta_{+-}=(2.232\pm 0.011)\times 10^{-3}$,
$\phi_{+-}=(43.51\pm 0.05)^{\circ}$. In the following calculations we neglect
the direct $CPV$ in kaons and assume $\varepsilon=\eta_{+-}$.

Despite the smallness of indirect CPV in kaons, itself it opens interesting
possibilities for flavour physics measurements Grossman:2011zk ; Bigi:1994aw ;
Kagan:2020vri . In the following section we consider few cases of charm hadron
decays that are particularly interesting for the $SU(3)_{f}$ probe.

##  III Strong phase difference between CF and DCS decays

Here we consider a set of two-body decays, where both Cabibbo-favourite and
doubly Cabibbo-suppressed amplitudes contribute to the final states. The
following analysis could be applied to the decays $D^{+}\to K^{0}_{S}\pi^{+}$,
$D_{s}^{+}\to K^{0}_{S}K^{+}$ and $\Lambda_{c}\to pK^{0}_{S}$ (the later is
not related to $SU(3)_{f}$ sum rules, but is interesting on its own).
Amplitudes of the final states with $K^{0}$ could be expressed then

$\displaystyle A_{f}$
$\displaystyle=\langle\pi\pi|H|\overline{K}{}^{0}\rangle+\sqrt{r_{f}}e^{i\delta}\langle\pi\pi|H|K^{0}\rangle,$
(18) $\displaystyle\overline{A_{f}}$
$\displaystyle=\langle\pi\pi|H|K^{0}\rangle+\sqrt{\overline{r_{f}}}e^{i\delta}\langle\pi\pi|H|\overline{K}{}^{0}\rangle$
(19)

where $\delta$, $r_{f}$ – strong phase difference and amplitude ratio for CF
and DCS amplitudes. In general case $r_{f}\neq\overline{r_{f}}$ and present
itself direct CPV in charm decays. Current world averages on $A_{CP}$ for
$D^{+}$ and $D_{s}^{+}$ mesons are Amhis:2019ckw :

$\displaystyle A_{CP}^{D_{s}\to K^{0}_{S}K^{+}}$ $\displaystyle=$
$\displaystyle(8\pm 26)\times 10^{-4},$ $\displaystyle A_{CP}^{D^{+}\to
K^{0}_{S}\pi^{+}}$ $\displaystyle=$ $\displaystyle(-41\pm 9)\times 10^{-4}.$
(20)

While non-zero effect was observed for $D^{+}$, it was noted Belle:2012ygt
that after subtraction $K^{0}$-$\overline{K}{}^{0}$-mixing contribution
asymmetry is consistent with zero. So we assume no CPV thereafter
$r_{f}=\overline{r_{f}}$.

Time-dependent decay rates could be obtained by substituting
$K^{0}$($\overline{K}{}^{0}$) decay amplitudes in (18, 19) with evolution
equations (14, 15).

$\displaystyle\mathcal{M}_{+}\equiv|A_{f}|^{2}=\mathcal{\overline{R}}$
$\displaystyle+$ $\displaystyle
r_{f}\mathcal{R}+\sqrt{r_{f}}\left(\cos{\delta}+2|\eta_{+-}|\sin{\delta}\sin{\phi_{+-}}\right)\times\Big{(}e^{-\Gamma_{S}t}-|\eta_{+-}|^{2}e^{-\Gamma_{L}t}\Big{)}$
$\displaystyle+$ $\displaystyle
2\sqrt{r_{f}}|\eta_{+-}|\Big{(}\sin{\delta+2|\eta_{+-}|\cos{\delta}\sin{\phi_{+-}}}\Big{)}e^{-\frac{1}{2}(\Gamma_{L}+\Gamma_{S})t}\sin{(\Delta
mt-\phi_{+-})},$
$\displaystyle\mathcal{M}_{-}\equiv|\overline{A_{f}}|^{2}=\mathcal{R}$
$\displaystyle+$ $\displaystyle
r_{f}\mathcal{\overline{R}}+\sqrt{r_{f}}\left(\cos{\delta}-2|\eta_{+-}|\sin{\delta}\sin{\phi_{+-}}\right)\times\Big{(}e^{-\Gamma_{S}t}-|\eta_{+-}|^{2}e^{-\Gamma_{L}t}\Big{)}$
$\displaystyle-$ $\displaystyle
2\sqrt{r_{f}}|\eta_{+-}|\Big{(}\sin{\delta}-2|\eta_{+-}|\cos{\delta}\sin{\phi_{+-}}\Big{)}e^{-\frac{1}{2}(\Gamma_{L}+\Gamma_{S})t}\sin{(\Delta
mt-\phi_{+-})}.$

These formulas demonstrate that the $K^{0}\to\pi^{+}\pi^{-}$ time-dependent
decay rates depend on the initial strong phase, moreover, both sine and cosine
of the strong phase enter the formula, therefore there are no trigonometrical
ambiguities in this measurement. The decay rates along with asymmetry are
illustrated in Fig. 2. One could see that the largest impact produced by the
strong phase on resulting asymmetry falls on big lifetimes
$\sim[6,14]\,\tau_{K_{S}}$.

Figure 2: Time-dependent decay rates for processes $D^{\pm}\to
K^{0}_{S}\pi^{\pm}$ (top) obtained with simulation. The bottom plots show the
resulting asymmetry and the difference between asymmetries with and without
DCS contribution. For this plot we used $\sqrt{r_{f}}=0.06$ and
$\delta^{0+}=-76^{\circ}$.

For the $D^{0}$-mesons situation gets complicates because of mixing. Time
evolution of flavour states is described by

$\displaystyle|D^{0}_{phys}(t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle
g_{+}(t^{\prime})\left|D^{0}\right\rangle-\left(\frac{q}{p}\right)_{\\!\\!D}g_{-}(t^{\prime})\left|\overline{D}{}^{0}\right\rangle,$
$\displaystyle|\overline{D}{}^{0}_{phys}(t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle
g_{+}(t^{\prime})\left|\overline{D}{}^{0}\right\rangle-\left(\frac{p}{q}\right)_{\\!\\!D}g_{-}(t^{\prime})\left|D^{0}\right\rangle,$
(23)

where $g_{\pm}=\frac{1}{2}\left(e^{-i\lambda_{2}t^{\prime}}\pm
e^{-i\lambda_{1}t^{\prime}}\right)$, $\lambda_{1,2}$ – are Hamiltonian
eigenvalues defined likewise eq. (12) and $t^{\prime}$ – $D^{0}$ decay time.
Here we assume no CPV in mixing, thus $|q/p|=1$. Amplitudes (18, 19) can be
rewritten as follows

$\displaystyle
A_{f}(t,t^{\prime})=a^{+}(t^{\prime})\langle\pi\pi|H|\overline{K}{}^{0}\rangle+b^{+}(t^{\prime})\langle\pi\pi|H|K^{0}\rangle,$
(24)
$\displaystyle\overline{A_{f}}(t,t^{\prime})=a^{-}(t^{\prime})\langle\pi\pi|H|K^{0}\rangle+b^{-}(t^{\prime})\langle\pi\pi|H|\overline{K}{}^{0}\rangle,$
(25)

where time-dependent coefficients are given by

$\displaystyle a^{+}(t^{\prime})$ $\displaystyle\equiv$
$\displaystyle\langle\overline{K}{}^{0}\pi^{0}|H|D^{0}_{phys}(t^{\prime})\rangle=A_{D^{0}}\left[g_{+}(t^{\prime})-\sqrt{r_{f}}e^{i(\delta+\phi)}g_{-}(t^{\prime})\right],$
$\displaystyle b^{+}(t^{\prime})$ $\displaystyle\equiv$ $\displaystyle\langle
K^{0}\pi^{0}|H|D^{0}_{phys}(t^{\prime})\rangle=A_{D^{0}}\left[\sqrt{r_{f}}e^{i(\delta-\phi)}g_{+}(t^{\prime})-g_{-}(t^{\prime})\right],$
(26)

and $a^{-}$, $b^{-}$ could be obtained from $a^{+}$, $b^{+}$ by substitution
$\phi\to-\phi$. Combined measurements of $D^{0}$-$\overline{D}{}^{0}$-mixing
yielded following values for mixing parameters Zyla:2020zbs :

$\displaystyle x$ $\displaystyle\equiv$ $\displaystyle\frac{\Delta
M}{\Gamma}=(0.43^{+0.10}_{-0.11})\%,$ (27) $\displaystyle y$
$\displaystyle\equiv$ $\displaystyle\frac{\Delta\Gamma}{2\Gamma}=(0.60\pm
0.06)\%,$ $\displaystyle\phi$ $\displaystyle\equiv$ $\displaystyle
Arg\left(\frac{q}{p}\right)_{\\!\\!D}=(0.08\pm 0.31)^{\circ}.$

In Kagan:2020vri  it was demonstrated that 2-dimensional distribution for the
$D^{0}\to K^{0}_{S}\pi^{0}$ decay is sensitive to a set of $CP$ observables.
However, it is unlikely that sufficient $D^{0}$-lifetime resolution will be
achieved in a neutral decay mode. Given the smallness of mixing parameters in
$D^{0}$-$\overline{D}{}^{0}$-system in the next section we consider a
possibility to integrate over $D^{0}$ lifetime and to use eq. (III, III) for
$D^{0}$/$\overline{D}{}^{0}$ as well.

##  IV Feasibility study

In this section we estimate potential precision of proposed measurement in
future experiments. It was shown in previous section that one would expect the
most sensitivity for strong phase could be achieved at big lifetimes –
$[6,14]\,\tau_{K_{S}}$. Based on this one could conclude that it is essential
that experiment should possess large tracking detector and/or produce soft
kaons. Also proper charged hadron identification is needed, since some of the
final states differ by $K/\pi$ interchange.

We consider the most promising experiments the project of future Super
$c$–$\tau$ factory Bondar:2013cja  and Belle II experiment Abe:2010gxa  that
is already taking data. Both experiments possess large drift chambers ($R\sim\
1$m) and produce relatively soft kaons, $\beta\gamma\sim 1..4$. Hadron
identification in Belle II provided with TOP in barrel part and ARICH in
endcaps and for the $c$–$\tau$ factory identification will be provided with
FARICH detector that covers almost full solid angle. Given the spatial
resolution of drift chambers $\sim 100\mu m$, kaon life time resolution could
be expected at the level of a few percent that is more than enough to perform
proposed measurement.

We perform feasibility study for the decay channels listed in Table 1. Future
Super $c$–$\tau$ factory is aiming to accumulate $10\mathrm{ab}^{-1}$ data
varying energies in c.m.s. from 3.097 GeV to 4.650 GeV. In particular
$3\mathrm{ab}^{-1}$ will be taken at $\psi(3770)$-resonance,
$1\mathrm{ab}^{-1}$ at $\psi(4160)$ and $1\mathrm{ab}^{-1}$ near
$\Lambda_{c}^{+}\Lambda_{c}^{-}$ threshold. For the Belle II experiment main
goal is $50\mathrm{ab}^{-1}$. To estimate potential yield of charm hadrons we
use $\sigma(ee\to c\bar{c})=1.1nb$ and fragmentation-fractions obtained in
Gladilin:2014tba ; Lisovyi:2015uqa . There is of course ambiguity due to event
selection criteria in each experiment and each particular channel. Here we
used conservative estimations for number of events, assuming only $30\%$ of
event will pass the selection for Belle II experiment and $70\%$ for
$c$–$\tau$ factory, since much cleaner environment is expected there. For the
$D^{0}$ studies we assumed $D^{*\pm}$ tagging in Belle II and semileptonic
tag-side decays for $\psi(3770)\to D^{0}\overline{D}{}^{0}$ case at Super
$c$–$\tau$ factory. Results are summarized in Table 1.

To confirm that there is no bias, we generate 100 MC samples of $40\times
10^{6}$ events, which correspond to $D^{+}\to K^{0}_{S}\pi^{+}$ decay, each
with a value of the angle $\delta$ in the $[-90^{\circ},90^{\circ}]$ interval
with a step of $10^{\circ}$. For the two-body decay modes with $K^{0}$ DCS/CF
amplitude ratios have not been measured yet, and could be approximated as

$r_{f}\equiv\left|\frac{\langle
K\pi|H|D\rangle}{\langle\overline{K}\pi|H|D\rangle}\right|^{2}\simeq\left|\frac{V_{cd}V_{us}^{*}}{V_{cs}V_{ud}^{*}}\right|^{2}\sim\mathcal{O}(\tan^{4}\theta_{c}).$
(28)

Analogous estimation for $D^{0}\to K^{+}\pi^{-}$ slightly differs from
experimental result, $r_{D}=(0.344\pm 0.002)\%$ Amhis:2019ckw .
$SU(3)_{f}$-breaking terms $K$ and $K^{\prime}$ introduced in (5, 6) aiming to
fix this small discrepancy on the amplitude level. Based on both experimental
and theoretical data it is reasonable to assume other DCS/CF ratios to be of
the same order $\sim\mathcal{O}(10^{-3})$. For this test we use
$\sqrt{r_{f}}=0.06$ which is very close to value measured for $D^{0}\to
K^{+}\pi^{-}$ decay. Each sample of MC contains time-dependent decay rates for
both particle and antiparticle (see example in Fig. 2). For each sample we
perform simultaneous unbinned maximum-likelihood fit for both time-dependent
decay rates. In the fit we consider events with $t/\tau_{K_{S}}>1$, since kaon
mixing does not contribute at low lifetimes. Fit results of one of the samples
presented in Fig. 3. Obtained results are in good agreement with generated
values of strong phase difference and amplitude ratio.

Table 1: Branching fraction and production yields |  | Estimated yield, $\times 10^{6}$ ($\delta$ uncertainty)
---|---|---
Channel |  Branching
fraction, $\%$ Zyla:2020zbs  |  Belle II ($50\ \mathrm{ab}^{-1}$) |  Super $c$–$\tau$
factory ($10\ \mathrm{ab}^{-1}$)
$D^{+}\to K^{0}_{S}\pi^{+}$ |  $1.56\pm 0.03$ |  40 ($5^{\circ}$) |  50 ($3^{\circ}$)
$D_{s}^{+}\to K^{0}_{S}K^{+}$ |  $1.46\pm 0.04$ |  20 ($7^{\circ}$) |  40 ($5^{\circ}$)
$\Lambda_{c}^{+}\to K^{0}_{S}p$ |  $1.59\pm 0.08$ |  15 ($8^{\circ}$) |  10 ($10^{\circ}$)
$D^{0}\to K^{0}_{S}\pi^{0}$ |  $1.23\pm 0.02$ |  30 ($6^{\circ}$) |  20 ($7^{\circ}$)
Figure 3: Results of feasibility study for the decay $D^{+}\to
K^{0}_{S}\pi^{+}$.

Since amplitude ratios $\sqrt{r_{f}}$ were not previously measured we perform
a scan for its values in the interval $[0.01,0.11]$ with step $2.5\times
10^{-3}$ simultaneously varying strong phase in the interval
$[-90^{\circ},90^{\circ}]$. Obtained uncertainties for 20 and 50 million
events, corresponding to $D^{0}\to K^{0}_{S}\pi^{0}$ and $D^{+}\to
K^{0}_{S}\pi^{+}$ decays shown in Fig. 4. As one could expect we observe
increasing sensitivity for strong phase with higher values of $\sqrt{r_{f}}$.
For the given amplitude ratio uncertainty in $\delta$ varies insignificantly
over the range $[-90^{\circ},90^{\circ}]$, which is certainly the advantage of
the method comparing it to usage of semileptonic $K^{0}$ decays
Pakhlov:2019wkn .

Measurements for the decay $D^{0}\to K^{0}_{S}\pi^{0}$ are of great importance
for the $SU(3)_{f}$ probe. Since achieving proper $D^{0}$ lifetime resolution
is hardly feasible we consider integration over $D^{0}$ lifetime –
$t^{\prime}$. For the purpose of the test MC generated distributions take into
account mixing effects, but fitting $p.d.f.$ are not. 1000 pseudo experiments
we performed and we found that amplitude ratio distribution turn out to be
shifted at about $1\sigma$ to the higher values. Such shift in general is
expected due to excess of “wrong”-flavour kaons arised from mixing. On the
other hand strong phase measurements still proved to be in good agreement with
generated values. We observed a $2^{\circ}$ bias in $\delta$, whereas the
statistical uncertainty is $6^{\circ}$.

Figure 4: Obtained uncertainties for the strong phase difference – $\delta$ in
$D^{+}$ decay (left) and in $D^{0}$ decay (right).

##  V Kaon regeneration

Proposed method based on the time-dependent study of $K^{0}\to\pi^{+}\pi^{-}$
decays, however beside flavour physics parameters some other effects could
contribute to decay rates. For example studies carries out in Ko:2010mk
showed that regeneration in the environment of today’s experiments could
induce a bias in the $A_{CP}$ measurement up to the level of $10^{-3}$.

In order to describe kaon propagation through matter the Hamiltonian in
Schrödinger equation (11) should be modified in the following way Good:1957zza

$i{\partial}_{t}{K^{0}(t)\choose\overline{K}{}^{0}(t)}=\left({{\mathbf{M}}-\frac{i}{2}{\mathbf{\Gamma}}}\right){K^{0}(t)\choose\overline{K}{}^{0}(t)}-\begin{pmatrix}\chi&0\\\
0&\overline{\chi}\end{pmatrix}{K^{0}(t)\choose\overline{K}{}^{0}(t)},$ (29)

where the second matrix describe nuclear scattering and coefficients defined
as

$\chi=-\frac{2\pi N}{m}f\ \ \mathrm{and}\ \ \overline{\chi}=-\frac{2\pi
N}{m}\overline{f},$ (30)

where $f(\overline{f})$ – are forward scattering amplitudes for
$K^{0}(\overline{K}{}^{0})$, $m$ – $K^{0}$ mass, $N=(\rho N_{A})/M$ – volume
density of the material, $N_{A}$ – Avogadro’s number, $\rho$ – mass density,
$M$ – molar mass. Strangeness conservation in strong interactions leads to
inequality forward scattering amplitudes $\Delta f\equiv f-\overline{f}\neq
0$. The evolution of the $K^{0}_{S}$($K^{0}_{L}$) beam could be expressed than

$\alpha_{S,L}=e^{-i\Sigma
t}\Big{[}\alpha_{S,L}^{0}\cos{\left(\frac{\Delta\lambda}{2}\sqrt{1+4r^{2}}t\right)}\pm
i\frac{\alpha_{S,L}^{0}\mp
2r\alpha_{L,S}^{0}}{\sqrt{1+4r^{2}}}\sin{\left(\frac{\Delta\lambda}{2}\sqrt{1+4r^{2}}t\right)}\Big{]},$
(31)

where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\lambda_{S}+\lambda_{L}+\chi+\overline{\chi}\right),$
$\displaystyle\Delta\lambda$ $\displaystyle=$
$\displaystyle\lambda_{S}-\lambda_{L},$ $\displaystyle\Delta\chi$
$\displaystyle=$ $\displaystyle\chi-\overline{\chi},$ $\displaystyle r$
$\displaystyle=$ $\displaystyle\frac{1}{2}\frac{\Delta\chi}{\Delta\lambda}.$
(32)

Regeneration parameter – $r$ is typically of the order of $10^{-2}$, so in the
following calculations we use the expansion for $\alpha_{S,L}$ to the lowest
order of $r$ (details could be found in Ref. Fetscher:1996fa ). It is
conventional to introduce the geometrical regeneration parameter:

$\zeta=r\left(1-e^{i\Delta\lambda\frac{Lm}{p}}\right),$ (33)

where $p$ – is kaon momentum and $L$ – regenerator thickness. Amplitudes (31)
could be expressed than in the form:

$\displaystyle\alpha_{S}(t)$ $\displaystyle=$ $\displaystyle
e^{\frac{1}{2}(\chi+\overline{\chi})t}e^{-i\lambda_{S}t}(\alpha_{S}^{0}+\zeta\alpha_{L}^{0}e^{-i\Delta\lambda
t}),$ $\displaystyle\alpha_{L}(t)$ $\displaystyle=$ $\displaystyle
e^{\frac{1}{2}(\chi+\overline{\chi})t}e^{-i\lambda_{L}t}(\alpha_{L}^{0}+\zeta\alpha_{S}^{0}).$
(34)

Applying equations (34) recursively for each passage through matter one could
account for kaon regeneration.

While accurate estimation of bias induced by regeneration should be performed
for each particular experiment, here we present an estimation based on typical
configurations. Since this study is mostly concerned with big kaon lifetimes,
we assume that neutral kaon have to pass through a beryllium beam pipe $(\sim
1\mathrm{mm})$ and a number of silicon layers of vertex detector. As a
reference Belle II configuration was used, where silicon vertex detector
consists of 6 layers ($L_{1,2}\simeq 50\mathrm{\mu m}$ and $L_{3-6}\simeq
300\mathrm{\mu m}$).

For this test we considered only leading regeneration contribution to CF decay
modes, since DCS/CF interference term is $\mathcal{O}(10^{-2})$ suppressed and
DCS term $\mathcal{O}(10^{-3})$ suppressed. We used the cross sections and
differences of forward scattering amplitudes obtained in Eberhard:1993nb
111Regeneration studies in CPLEAR experiment CPLEAR:1997fwa  showed good
agreement between optical model predictions and experimental results. Using MC
simulation we found that for $1\mathrm{GeV/c}$ kaons bias in the strong phase
measurement is under $4^{\circ}$. Obtained value is comparable with potential
statistical uncertainty, however regeneration could be the main source of
systematic uncertainty and for each particular environment studies are
required.

##  VI Summary

In this paper we presented a method to measure strong phase differences in
charm hadron decay with $K^{0}$-meson in the final state. It was shown that
CPV in $K^{0}$-$\overline{K}{}^{0}$ mixing allows us to disentangle initial
combination
$a\left|K^{0}(t)\right\rangle+b\left|\overline{K}{}^{0}(t)\right\rangle$ that
arises in the presence of CF and DCS decays. In order to perform such
measurement the experiment should satisfy following requirements: large
tracking detector that allows reconstruct $K^{0}_{S}$ decays even after
$10\tau_{K_{S}}$, sufficient statistics of charm decays –
$\mathcal{O}(10^{6}..10^{7})$ and proper charged hadron identification. The
Belle II experiment and future Super $c$–$\tau$ factory are good candidates
for such measurement. While LHC experiments have huge data samples of charmed
mesons, some of the features of the detectors significantly reduce the
possibilities of such measurements: LHCb has too short tracker, while CMS and
Atlas have no particle identification.

To estimate potential precision of the method feasibility study was performed.
Expected number of events was calculated for each particular channel for both
experiments. Proposed measurements proved to be unbiased and free of
trigonometrical ambiguity. Obtained results for statistical uncertainty
(assuming $r_{f}\sim\mathcal{O}(10^{-3})$) are comparable with current
theoretical uncertainties and uncertainties that could be obtained with
semileptonic kaon decays. We also presented an estimation of regeneration
contributions to proposed measurement.

##  VII Acknowledgments

The work of P. Pakhlov was conducted within the framework of the Basic
Research Program at the National Research University Higher School of
Economics (HSE). The reported study of V. Popov was funded by RFBR, project
number 19-32-90104.

## References

  * (1) Y. Grossman and D. J. Robinson, “SU(3) Sum Rules for Charm Decay,” JHEP 1304, 067 (2013)
  * (2) R. Aaij et al. [LHCb], “Observation of CP Violation in Charm Decays,” Phys. Rev. Lett. 122, no.21, 211803 (2019)
  * (3) M. Golden and B. Grinstein, “Enhanced CP Violations in Hadronic Charm Decays,” Phys. Lett. B 222, 501-506 (1989)
  * (4) G. Hiller, M. Jung and S. Schacht, “SU(3)-flavor anatomy of nonleptonic charm decays,” Phys. Rev. D 87, no.1, 014024 (2013)
  * (5) J. Brod, Y. Grossman, A. L. Kagan and J. Zupan, “A Consistent Picture for Large Penguins in D -$>$ pi+ pi-, K+ K-,” JHEP 10, 161 (2012)
  * (6) E. Franco, S. Mishima and L. Silvestrini, “The Standard Model confronts CP violation in $D^{0}\to\pi^{+}\pi^{-}$ and $D^{0}\to K^{+}K^{-}$,” JHEP 05, 140 (2012)
  * (7) F. Buccella, A. Paul and P. Santorelli, “$SU(3)_{F}$ breaking through final state interactions and $CP$ asymmetries in $D\to PP$ decays,” Phys. Rev. D 99, no.11, 113001 (2019)
  * (8) P. Pakhlov and V. Popov, “Measurement of $D^{0}-\bar{D}^{0}$ mixing parameters using semileptonic decays of neutral kaon,” JHEP 02, 160 (2020)
  * (9) P.A. Zyla et al. [Particle Data Group], “Review of Particle Physics,” PTEP 2020, no.8, 083C01 (2020)
  * (10) Y. Grossman and Y. Nir, “CP Violation in $\tau^{\pm}\to\pi^{\pm}K_{S}\nu$ and $D^{\pm}\to\pi^{\pm}K_{S}$: The Importance of $K_{S}-K_{L}$ Interference,” JHEP 04, 002 (2012)
  * (11) I. I. Y. Bigi and H. Yamamoto, “Interference between Cabibbo allowed and doubly forbidden transitions in D —$>$ K(S), K(L) + pi’s decays,” Phys. Lett. B 349, 363-366 (1995)
  * (12) A. L. Kagan and L. Silvestrini, “Dispersive and absorptive $CP$ violation in $D^{0}-\overline{D^{0}}$ mixing,” Phys. Rev. D 103, no.5, 053008 (2021)
  * (13) Y. S. Amhis et al. [HFLAV], “Averages of b-hadron, c-hadron, and $\tau$-lepton properties as of 2018,” Eur. Phys. J. C 81, no.3, 226 (2021)
  * (14) B. R. Ko et al. [Belle], “Evidence for CP Violation in the Decay $D^{+}\rightarrow K^{0}_{S}\pi^{+}$,” Phys. Rev. Lett. 109, 021601 (2012) [erratum: Phys. Rev. Lett. 109, 119903 (2012)]
  * (15) T. Abe et al. [Belle-II], “Belle II Technical Design Report,” [arXiv:1011.0352 [physics.ins-det]].
  * (16) E. Kou et al. [Belle-II], “The Belle II Physics Book,” PTEP 2019, no.12, 123C01 (2019) [erratum: PTEP 2020, no.2, 029201 (2020)]
  * (17) A. E. Bondar et al. [Charm-Tau Factory Collaboration], “Project of a Super Charm-Tau factory at the Budker Institute of Nuclear Physics in Novosibirsk,” Phys. Atom. Nucl. 76, 1072 (2013).
  * (18) M. Lisovyi, A. Verbytskyi and O. Zenaiev, “Combined analysis of charm-quark fragmentation-fraction measurements,” Eur. Phys. J. C 76, no.7, 397 (2016)
  * (19) L. Gladilin, “Fragmentation fractions of $c$ and $b$ quarks into charmed hadrons at LEP,” Eur. Phys. J. C 75, no.1, 19 (2015)
  * (20) B. R. Ko, E. Won, B. Golob and P. Pakhlov, “Effect of nuclear interactions of neutral kaons on CP asymmetry measurements,” Phys. Rev. D 84, 111501 (2011)
  * (21) M. L. Good, “Relation between Scattering and Absorption in the Pais-Piccioni Phenomenon,” Phys. Rev. 106, 591-595 (1957)
  * (22) W. Fetscher, P. Kokkas, P. Pavlopoulos, T. Schietinger and T. Ruf, “Regeneration of arbitrary coherent neutral kaon states: A new method for measuring the K0 anti-K0 forward scattering amplitude,” Z. Phys. C 72, 543-547 (1996)
  * (23) P. H. Eberhard and F. Uchiyama, “Neutral kaon regeneration probabilities at asymmetric phi factory energies,” Nucl. Instrum. Meth. A 350, 144-149 (1994)
  * (24) A. Angelopoulos et al. [CPLEAR], “Measurement of the neutral kaon regeneration amplitude in carbon at momenta below 1-GeV/c,” Phys. Lett. B 413, 422-430 (1997)