# Dressed Majorana fermion in a hybrid nanowire

Guo-Jian Qiao  Beijing Computational Science Research Center, Beijing 100193,
China  Xin Yue  Beijing Computational Science Research Center, Beijing 100193,
China  C. P. Sun<EMAIL_ADDRESS>Beijing Computational Science Research
Center, Beijing 100193, China  Graduate School of China Academy of Engineering
Physics, Beijing 100193, China  School of Physics, Peking University, Beijing
100871, China

###### Abstract

The low-energy theory of hybrid nanowire systems fails to define Majorana
fermion (MF) in the strong tunneling and magnetic field strength. To address
this limitation, we propose a holistic approach to define MF in which the
quasi-excitation in nanowire and superconductor constitutes together its own
“antiparticles”. This definition is general, beyond the constraint presented
in the low-energy theory. It reveals that the Majorana phase depends not only
on the chemical potential and Zeeman energy in nanowire but also on those of
superconductor, and that the mismatch of chemical potential leads not to
observe MF. Such a broader perspective provides more specific experimental
guidance under various conditions.

Introduction.—The observation of Majorana fermion (MF) in condensed matter
systems is a prerequisite for realizing topological quantum computation
(Kitaev, 2001; Sau  _et al._ , 2010; Alicea, 2012; Lutchyn  _et al._ , 2018;
Prada  _et al._ , 2020). As one of the promising platforms, the semiconductor
nanowire with appreciable spin-orbit coupling coupled to an s-wave
superconductor (SC), as considered to display MF (Lutchyn  _et al._ , 2010;
Oreg  _et al._ , 2010; Mourik  _et al._ , 2012), mainly based on a low-energy
effective theory (Alicea, 2012; Stanescu  _et al._ , 2011; Stanescu and Das
Sarma, 2017). Such a low-energy effective theory predicts that the topological
phase of MF can be implemented in experiments by only adjusting the chemical
potential and applied magnetic field in the nanowire, directly independent on
the proximity superconducting materials.

The above low-energy effective theory is derived by eliminating the virtual
process of exchange of electron in nanowire between Bogoliubov quasi-particles
in SC, and a constant s-wave pairing is induced in the nanowire. But the
validity of the low-energy theory is limited due to the following three
approximations (Stanescu  _et al._ , 2011; Stanescu and Das Sarma, 2017; sup,
):

(i) the s-wave SC is assumed at a periodic boundary to ignore the boundary
effect, which indicates that the edge state of MF is only localized at two
ends of the nanowire;

(ii) the energy-independent coupling spectral between the nanowire and SC in
wide-band limit results in the induced constant pairing term;

(iii) only low-energy effects of superconducting self-energy correction are
considered, which limits the effective theory to the weak tunneling strength
between the nanowire and SC.

What is more, it has also been pointed out that the effective Kitaev model
cannot be reduced when the Zeeman splitting and tunneling strength are
comparable to the superconducting gap due to the failure of perturbation
method (Qiao  _et al._ , 2022), which leaves unresolved the question whether
MF can be defined in the strong magnetic field and coupling region.

Figure 1: The semiconductor nanowire is partially or fully covered by the
superconducting shell, which can be characterized by the one-dimensional
lattice model of nanowire coupled to an s-wave superconductor. And the
spectrum of its differential conductance is obtained by connecting electron
leads at the ends of the nanowire.

However, in the mainstream hybrid system used to observe MF experimentally
such as InSb/InAs nanowire partially or fully covered by aluminum film, the
coupling strength between nanowire and SC has reached the intermediate or even
strong coupling schemes under the strong magnetic field (Das  _et al._ , 2012;
Valentini  _et al._ , 2021; Zhang  _et al._ , 2018; Vaitiekenas  _et al._ ,
2020; Lutchyn  _et al._ , 2018). Besides, some artificially experimental
signature of MF fitting a phase diagram given by the above low-energy theory
have caused great controversy and were even retracted so that the validity of
the low-energy theory should be doubted (Zhang  _et al._ , 2018; Gazibegovic
_et al._ , 2017; Vaitiekenas  _et al._ , 2020; He  _et al._ , 2017).
Therefore, it is necessary to address the following issues:

(i) how to generally define MF in a hybrid system at any magnetic field and
tunneling strength ?

(ii) what is universal Majorana topological phase like? Does it also depend on
the relevant parameters of SC ?

(iii) under what conditions, the topological phase of MF obtained by the low-
energy effective Hamiltonian is valid ?

In this paper, by a minimal hybrid model (Lutchyn  _et al._ , 2018) that a
one-dimensional semiconductor nanowire couples to a one-dimensional s-wave SC
[Fig. 1], we propose a holistic perspective to define the dressed Majorana
fermion (DMF), namely, DMF is composed of the electron and hole excitations in
nanowire dressed by quasi-excitations in SC. Here, we treat the excitations in
nanowire and SC equally, rather than directly eliminating quasi-excitation in
SC, as has been done in the low-energy effective theory. Therefore, such a
definition of DMF holds for any tunneling and magnetic field strength.

The edge state of DMF in open boundary is localized not only at two ends of
nanowire but also of SC, and its distribution in nanowire and SC is determined
by the tunneling strength. As a signature of DMF, we show that the $2e^{2}/h$
zero-bias peak in the differential conductance spectrum is also obtained by
connecting electron leads to nanowire.

We demonstrate that the universal Majorana phase depends not only on the
chemical potential and Zeeman energy in nanowire but also in SC. The low-
energy phase diagram given by the effective theory [it depends only on
chemical potential and magnetic field in nanowire] is only an approximate
result of the universal Majorana phase in the weak tunneling scheme.

Lastly, the universal Majorana phase is applied to the experimental hybrid
InSb/InAs-Al system, and it is suggested that the mismatch of chemical
potentials between nanowire and SC hinders the observation of MF.

Theoretical Model.—The semiconductor nanowire with the Rashba spin-orbit
coupling $\alpha_{w}$ is described by $N$-sites lattice model as (Stanescu
_et al._ , 2011; Stoudenmire  _et al._ , 2011; Chevallier and Klinovaja, 2016)

$\displaystyle H_{w}$
$\displaystyle=\sum_{n=1}^{N}d_{n}^{\dagger}\,[(t_{w}-\mu_{w})\sigma_{0}+h_{w}\sigma_{x}]\,d_{n}$
$\displaystyle-\sum_{n=1}^{N-1}d_{n}^{\dagger}\,(\frac{t_{w}}{2}\sigma_{0}+\frac{i\alpha_{w}}{2}\sigma_{y})\,d_{n+1}+\mathrm{H.c.},$
(1)

where $d_{n}:=[d_{n\uparrow},d_{n\downarrow}]^{T}$ is a annihilation operator
of the electron in site-$n$ with spin $\uparrow,\downarrow$, and
$\sigma_{0,x,y}$ are the $2\times 2$ identity and Pauli matrices. Here,
$\mu_{w}$ is the chemical potential, $t_{w}$ is the hopping strength between
nearest sites in the nanowire, and $h_{w}$ is the Zeeman splitting caused by
an external magnetic aligned with the nanowire axis.

The semiconductor nanowire is in contact with an s-wave SC or wrapped by
superconducting shell (such as Al (Das  _et al._ , 2012; Zhang  _et al._ ,
2018) and Pb (Cao  _et al._ , 2022; Jiang  _et al._ , 2022) shell), and such a
SC is also characterized by the one-dimensional lattice model: (Reeg  _et al._
, 2017)

$\displaystyle H_{s}$
$\displaystyle=\sum_{n=1,\sigma}^{N}(t_{s}-\mu_{s})c_{n,\sigma}^{\dagger}c_{n,\sigma}-\frac{t_{s}}{2}[c_{n,\sigma}^{\dagger}c_{n+1,\sigma}+\mathrm{H.c.}]$
$\displaystyle+\sum_{n=1}^{N}[\Delta_{s}c_{n,\uparrow}^{\dagger}c_{n,\downarrow}^{\dagger}+h_{s}c_{n\uparrow}^{\dagger}c_{n\downarrow}+\mathrm{H.c.]}.$
(2)

Similarly, $\mu_{s}$ is the chemical potential, $t_{s}$ is hopping strength,
$h_{s}$ is the Zeeman energy (Clogston, 1962; Maki and Tsuneto, 1964; Reeg and
Maslov, 2017; Vaitiekėnas  _et al._ , 2018) and $\Delta_{s}$ is a pairing
strength for the s-wave SC.

And the interaction between nanowire and SC is via single electron tunneling:
$H_{t}=-\sum_{n,m,\sigma}T_{nm}[d_{n\sigma}c_{m\sigma}^{\dagger}+c_{m\sigma}d_{n\sigma}^{\dagger}],$
with the tunneling strength $T_{nm}$ between the site-$n$ of the nanowire and
site-$m$ in SC. Below, we consider that the tunneling strength for the same
sites $n=m$ is equal but ignored between different sites $n\neq m$, i.e.,
$T_{nm}\equiv\delta_{mn}T$.

Notice that the Hamiltonian of the hybrid nanowire can be diagonalized as
$H=H_{w}+H_{s}+H_{t}=1/2\sum_{E}E\,\eta_{E}^{\dagger}\eta_{E}$, and the
corresponding quasi-particle operator is

$\eta_{E}=(\mathbf{u}_{E}^{w})^{\dagger}\cdot\mathbf{d}+(\mathbf{v}_{E}^{w})^{\dagger}\cdot(\mathbf{d}^{\dagger})^{T}+(\mathbf{u}_{E}^{s})^{\dagger}\cdot\mathbf{c}+(\mathbf{v}_{E}^{s})^{\dagger}\cdot(\mathbf{c}^{\dagger})^{T},$
(3)

where $\mathbf{d}:=[d_{1},\ldots,d_{N}]^{T}$ and
$\mathbf{c}:=[c_{1},\ldots,c_{N}]^{T}$ with
$c_{n}:=[c_{n\uparrow},c_{n\downarrow}]^{T}$. And the eigen-wave function is
$\Psi_{E}=[\mathbf{u}_{E},\mathbf{v}_{E}]^{T}$ with
$\mathbf{u}_{E}=[\mathbf{u}_{E}^{w},\mathbf{u}_{E}^{s}]^{T}$ and
$\mathbf{v}_{E}=[\mathbf{v}_{E}^{w},\mathbf{v}_{E}^{s}]^{T}$, where
$\mathbf{o}_{E}^{\alpha}=[\mathrm{o}_{E,1}^{\alpha},\ldots,\mathrm{o}_{E,N}^{\alpha}]^{T}$
with
$\mathrm{o}_{E,n}^{\alpha}=[\mathrm{o}_{E,n\uparrow}^{\alpha},\mathrm{o}_{E,n\downarrow}^{\alpha}]^{T}$
respctively represent the wave function of electrons
$\mathrm{o}\equiv\mathrm{u}$ and holes $\mathrm{o}\equiv\mathrm{v}$ in the
nanowire and SC for $\alpha=w,s$.

Dressed MF.—We define the dressed Majorana fermion as

$\eta_{E}=\eta_{E}^{\dagger}\Longrightarrow\mathbf{u}_{E}^{w}=(\mathbf{v}_{E}^{w})^{*},\;\mathbf{u}_{E}^{s}=(\mathbf{v}_{E}^{s})^{*},$
(4)

where the quasi-excitation constitutes its own “antiparticles”. Further, the
particle-hole symmetry requires $(\mathbf{u}_{-E}^{w})^{*}=\mathbf{v}_{E}^{w}$
and $(\mathbf{u}_{-E}^{s})^{*}=\mathbf{v}_{E}^{s}$. Thus it is proved that
only the zero-energy quasi-particle is Majorana fermion (sup, ). Then, we need
to give the conditions for the existence of such a zero-energy wave function.

For the zero-enrgy wave function $E=0$ ($\mathbf{u}=\mathbf{v}^{*}$, where
subscript “$0$” is omitted), we decompose $\mathbf{u}$ into real and imaginary
parts $\mathbf{u}=\mathbf{u}^{(r)}+i\mathbf{u}^{(i)}$. It is easy to verify
that the electron wave function’s real and imaginary parts in nanowire satisfy

$[(\mathbf{H}_{s}+\lambda\mathbf{P}_{s})\cdot\mathbf{H}_{w}-\mathbf{T}^{2}]\cdot\mathbf{u}_{w}^{(\lambda)}=0,$
(5)

and the electron wave function in SC is
$\mathbf{u}_{s}^{(\lambda)}=-\mathbf{T}^{-1}\cdot\mathbf{H}_{w}\cdot\mathbf{u}_{w}^{(\lambda)}$.
Here, $\lambda=\pm 1$ corresponds to the real and imaginary parts of
$\mathbf{u}_{\alpha}^{(\lambda)}$, that is
$\mathbf{u}_{\alpha}^{(1)}\equiv\mathbf{u}_{\alpha}^{(r)}$ and
$\mathbf{u}_{\alpha}^{(-1)}\equiv\mathbf{u}_{\alpha}^{(i)}$ with $\alpha=w,s$.
And the $2N\times 2N$ Hamiltonian matrices of the nanowire $\mathbf{H}_{w}$,
tunneling interaction $\mathbf{T}$ and SC ($\mathbf{H}_{s}$ and
$\mathbf{P}_{s}$) are given in supplemental materials (sup, ).

We assume that
$\mathrm{u}_{w,n}^{(r)}\equiv\xi^{n}[a_{\uparrow},a_{\downarrow}]^{T},1\leq
n\leq N$, where $\xi$ is a complex number, and $a_{\uparrow},a_{\downarrow}$
are undetermined coefficients. By substituting $\mathrm{u}_{w,n}^{(r)}$ into
Eq. (5), it is proved that $\xi$ is root of $8$-order equation of one variable
$f_{\lambda}(\xi)=0$ with $\lambda=1$, and $a_{\uparrow},a_{\downarrow}$ are
determined for an given $\xi$ (sup, ). Thus such a electron wave function in
nanowire becomes
$\mathrm{u}_{w,n}^{(r)}=\sum_{j=1}^{8}\alpha_{j}\xi_{j}^{n}[a_{j\uparrow},a_{j\downarrow}]^{T},$
where $\alpha_{j}$ are arbitrary superposition coefficients.

The boundary and normalization conditions of the zero-energy wave function
$\mathbf{u}^{(r)}$ require that $f_{1}(\xi)=0$ must have three roots less than
$1$ as the total site number approaches infinity $N\rightarrow\infty$, namely

$f_{1}(\xi_{j})=0,\exists\>\xi_{j},j=1,2,3,\>|\xi_{j}|<1$ (6)

At this time, $\mathbf{u}_{w}^{(r)}$ is uniquely determined, and the wave
function in SC is obtained by
$\mathbf{u}_{s}^{(r)}=-\mathbf{T}^{-1}\cdot\mathbf{H}_{w}\cdot\mathbf{u}_{w}^{(r)}$.

Similarly, the imagine part of the wave function $\mathbf{u}_{w}^{(i)}$ is
obtained as
$\mathrm{u}_{w,n}^{(i)}=\sum_{j=1}^{3}\beta_{j}\xi_{j}^{-N+n+1}[b_{j\uparrow},b_{j\downarrow}]^{T}$
for the given $\beta_{j}$ and $b_{j\sigma}$ with $\sigma=\pm 1$, and
$\mathbf{u}_{s}^{(i)}=-\mathbf{T}^{-1}\cdot\mathbf{H}_{w}\cdot\mathbf{u}_{w}^{(i)}$.
Here, we have utilized the function relation
$f_{\lambda=1}(\xi)=f_{\lambda=-1}(1/\xi)=0$ with $\xi\neq 0$ (sup, ).

Known $\mathbf{u}^{(r)}$ and $\mathbf{u}^{(i)}$, DMF is rewritten as

$\displaystyle\gamma_{1}$
$\displaystyle=\sum_{n\sigma}\mathrm{u}_{w,n\sigma}^{(r)}(d_{n\sigma}+d_{n\sigma}^{\dagger})+\mathrm{u}_{s,n\sigma}^{(r)}(c_{n\sigma}+c_{n\sigma}^{\dagger}),$
(7) $\displaystyle\gamma_{N}$
$\displaystyle=\sum_{n\sigma}\mathrm{u}_{w,n\sigma}^{(i)}(i[d_{n\sigma}^{\dagger}-d_{n\sigma}])+\mathrm{u}_{s,n\sigma}^{(i)}(i[c_{n\sigma}^{\dagger}-c_{n\sigma}]).$

It is seen from Eqs. (6, 7) that $\gamma_{1}$ and $\gamma_{N}$ correspond to
the localized edge state at site-$1$ and site-$N$ of the nanowire and SC. Such
a MF is a compsite of electron and hole exication in nanowire dressed by the
quasi-exication of SC [the last two terms of (7)], and it is different from MF
defined by the low-energy effective theory, where the quasi-excitation in SC
has been eliminated.

Apparently, there is another case where the real and imaginary part of zero-
energy wave function are respetively localized near site-$N$ and site-$1$ of
the nanowire and SC [$|\xi_{j}|>1$ in Eq. (6)]. In short, a zero-energy wave
function exists only when the function $f_{1}(\xi)$ has three complex roots,
and their magnitudes are greater than or less than $1$. Thus, the parameter
range in the presence of the edge state of DMF is

$\prod_{k}[(h_{w}+Z_{k}h_{s})^{2}-(\epsilon_{k}^{w}-\mu_{w}-Z_{k}(\epsilon_{k}^{s}-\mu_{s}))^{2}-(Z_{k}\Delta_{s})^{2}]<0,$
(8)

where $\epsilon_{k}^{s}:=t_{\alpha}(1-\cos k)$, and $k$ only takes $0$ and
$\pi$. The correction factor $Z_{k}:=T^{2}/(E_{k,+}^{s}E_{k,-}^{s})$ embodies
the dressed effect of SC for the nanowire, where the excitation energy of
quasi-particle in SC is included as
$E_{k,\pm}^{s}=\sqrt{[\epsilon_{k}^{s}-\mu_{s}]^{2}+\Delta_{s}^{2}}\pm h_{s}$
(Alicea, 2012). Such a edge state of DMF depends not only on the chemical
potential $\mu_{w}$ and Zeeman energy $h_{w}$ of the nanowire, but also on the
chemical potential $\mu_{s}$ and Zeeman energy $h_{s}$ in SC. And it can be
proved that the conditions for the existence of the Majorana edge state (8) is
consistent with the topological phase of DMF (Majorana phase) given by the
topological invariant
$\text{sgn}(\mathrm{Pf}[\mathscr{H}_{\gamma}(0)])\text{sgn}(\mathrm{Pf}[\mathscr{H}_{\gamma}(\pi)])<0$
(sup, ).

In this Majorana phase, the energy spectrum of hybrid system indeed displays
the zero-energy modes of DMF, and there is a gap between the zero and
excitation energy [Fig. 2(a)]. The edge state of DMF is not only localized at
two ends in the nanowire, but also in SC [Fig. 2(c, d)]. And the proportion of
zero-energy wave function distributed in the nanowire and SC
$\mathrm{P}_{\mathrm{SC}}/\mathrm{P}_{\mathrm{NW}}$ increases with the enhance
of tunneling strength $T$ , where the zero-energy wave function of the
nanowire and SC are measured by
$\mathrm{P}_{\mathrm{NW}}=\sum_{n\sigma}|\mathrm{u}_{w,n\sigma}|^{2}+|\mathrm{v}_{w,n\sigma}|^{2}$
and
$\mathrm{P}_{\mathrm{SC}}=\sum_{n\sigma}|\mathrm{u}_{s,n\sigma}|^{2}+|\mathrm{v}_{s,n\sigma}|^{2}$
with normalization condition
$\mathrm{P}_{\mathrm{SC}}+\mathrm{P}_{\mathrm{NW}}=1$ [Fig. 2(b)]. This shows
that the more zero-energy wave function penetrates SC from the nanowire as the
tunneling strength increases. However, the low-energy effective theory has
ignored the zero-energy wave function in SC since the superconductor is
considered as a periodic boundary. In the intermediate, especially strong
tunneling schemes $T>\Delta_{s}$, the zero-energy wave function in SC is
comparable to that in nanowire [$T\sim
2\Delta_{s},\>\mathrm{P}_{\mathrm{SC}}/\mathrm{P}_{\mathrm{NW}}\sim 30\%$],
and the boundary effect of SC cannot be ignored.

Figure 2: (a) The zero-energy and excitation-energy modes in the hybrid
system. (b) The zero-energy wave function
$\mathrm{P}_{\mathrm{SC}},\mathrm{P}_{\mathrm{NW}}$ in nanowire and SC, and
their proportion $\mathrm{P}_{\mathrm{SC}}/\mathrm{P}_{\mathrm{NW}}$ change as
the tunneling strength enhances. (c, d) The distribution of the localized
zero-mode wave functions at each site-$n$ in nanowire and SC. (e) The zero
bias peak in differential spectrum by connecting electron leads at the ends of
the nanowire. The parameters are set as
$N=600,t_{w}=12\Delta_{s},t_{s}=10\Delta_{s},\mu_{w}=0,\mu_{s}=4\Delta_{s},h_{w}=1.5\Delta_{s},h_{s}=0,\alpha_{w}=T=1.5\Delta_{s}.$
(f) The general Majorana phase (light purple region) deviates from the phase
diagram determined by the low-energy theory (red dashed line) in
$\mu_{w}-h_{w}$ space with the correction factor $Z_{0}=0.2$ and
$h_{s}=0.1h_{w}$.

It is known that MF determined by the low-energy theory will result in the
zero-bias peak (ZBP) with the height $2e^{2}/h$ at zero temperature (Qiao  _et
al._ , 2022; Li  _et al._ , 2014; Flensberg, 2010). For DMF in the hybrid
nanowire, the $2e^{2}/h$ ZBP as a signature of DMF also appears by connecting
electron leads at the ends of the nanowire [see Fig. (1) and Fig. 2(e)] (sup,
).

Typically, the bandwidth of the nanowire and SC as the maximum energy scale is
much larger than their Zeeman energy, chemical potential and the
superconducting gap $t_{\alpha}\gg\mu_{\alpha},h_{\alpha},\Delta_{s}$ with
$\alpha=w,s$. Then the universal Majorana phase (8) becomes

$(h_{w}+Z_{0}h_{s})>\sqrt{(\mu_{w}-Z_{0}\mu_{s})^{2}+(Z_{0}\Delta_{s})^{2}}.$
(9)

When the dressed effect of SC for the nanowire (the correction factor $Z_{0}$)
is small so that $Z_{0}h_{s},Z_{0}\mu_{s}\ll 1$, the Majorana phase returns to
the low-energy phase diagram $h_{w}>\sqrt{\mu_{w}^{2}+\Delta^{2}}$ with
$\Delta:=Z_{0}\Delta_{s}$ (Alicea, 2012; Lutchyn  _et al._ , 2018). In this
regard, the correction factor can be regraded as the ratio of the induced
pairing strength $\Delta$ to the superconducting gap
$Z_{0}=\Delta/\Delta_{s}$, which characterizes the different tunneling
schemes. By increasing the tunneling strength between nanowire and SC, the
correction factor gradually achieves a stronger tunneling scheme
$Z_{0}\rightarrow 1$ so that the Zeeman energy
$\tilde{h}_{w}:=h_{w}+Z_{0}h_{s}$ and chemical potential
$\tilde{\mu}_{w}:=\mu_{w}-Z_{0}\mu_{s}$ in the nanowire is significantly
modified. This results in the Majorana phase deviating significantly from the
previous low-energy phase diagram, especially for the more considerable
superconducting chemical potentials and Zeeman energy [Fig. 2(f)]. Therefore,
the Majorana phase (9) is more general, the fragile low-energy topological
phase is just an approximate result of which in weak tunneling scheme
$Z_{0}\ll 1$.

It is worth noting that the low-energy phase diagram depends only on chemical
potential $\mu_{w}$ and magnetic field $h_{w}$ in nanowire, which is
considered as the solid theoretical basis to observe MF in experiments so far.
However, the universal Majorana phase also depends significantly on the
chemical potential $\mu_{s}$ and magnetic field $h_{s}$ in SC [see Eq. (9)],
and it does have an impact on the observation of MF in an actual hybrid
nanowire system.

Dependence of Majorana phase on the chemical potential and magnetic field in
SC—For the hybrid InSb/InAs-Al system, the superconducting gap of aluminum
film is $\Delta_{s}=0.34\,\mathrm{meV}$, and the Zeeman splitting in SC is an
order of magnitude smaller than that in the nanowire due to the difference in
Landé factor $h_{s}\sim 0.05-0.15h_{w}$ (Cao  _et al._ , 2022; Reeg and
Maslov, 2017; Vaitiekėnas  _et al._ , 2018; Lutchyn  _et al._ , 2018). From
the induced energy gap observed in experiments $\Delta=0.2\,\mathrm{meV}$
(Zhang  _et al._ , 2018), the correction factor can be estimated as
$Z_{0}=\Delta/\Delta_{s}\simeq 0.6$. Therefore the magnetic field required for
the Majorana phase by (9) is obtained as

$B>\frac{2\sqrt{(\mu_{w}-0.6\mu_{s})^{2}+(0.6\Delta_{s})^{2}}}{\mu_{\mathrm{B}}(\mathrm{g}_{w}+0.6\mathrm{g}_{s})}.$
(10)

with the Bohr magneton $\mu_{\mathrm{B}}$ and the Landé factor
$\mathrm{g}_{\alpha},\alpha=w,s$ of nanowire and SC. Since the critical
magnetic field of the aluminum film is $2\,\mathrm{T}$ (Cao  _et al._ , 2022),
the maximum difference in chemical potential between the nanowire and SC
cannot be greater than $|\mu_{w}-0.6\mu_{s}|<2.4/0.9\,\mathrm{meV}$ for the
InSb\InAs nanowire, where the Landé factors are taken as
$\mathrm{g}_{\mathrm{InSb}}=40,\mathrm{g}_{\mathrm{InAs}}=15$ and
$\mathrm{g}_{s}=2$ (Cao  _et al._ , 2022; Lutchyn  _et al._ , 2018). For
example, if the superconducting chemical potential is $1\mathrm{eV}$
(Mikkelsen  _et al._ , 2018), the chemical potential of the InSb\InAs nanowire
must be adjusted to $597/599\lesssim\mu_{w}\apprle 602/600\>\mathrm{meV}$.
Only within such a small range of the chemical potential, MF is possibly
observed. The smaller the external magnetic field, the smaller the above
window of chemical potential of nanowire for the Majorana phase. Therefore,
the mismatch of the chemical potentials in nanowire and SC will prevent the
Majorana topological phase. This may understand why the signature of zero-bias
peak in the current hybrid nanowire systems is non-topological (Valentini  _et
al._ , 2021; Zhang  _et al._ , 2018; Wang  _et al._ , 2022).

Based on the above facts, there are some remarks for observing MF in
experiments. (i) The magnetic field should be modest strength. An excessively
large magnetic field will inhibit the superconducting gap (Liu  _et al._ ,
2017; Mohapatra  _et al._ , 2019; Pan and Das Sarma, 2020; Pan  _et al._ ,
2024), while a smaller magnetic field will require higher control accuracy of
the chemical potential in nanowire and SC. (ii) The mismatch of chemical
potential highlighted above requires that the chemical potential in the
nanowire is adjusted in an extensive range to achieve the Majorana phase,
which may produce the current between the nanowire and SC due to the drastic
change of their potential difference. This is not conducive to observe MF.
(iii) The instability of (especially the larger) superconducting chemical
potential as other experimental parameters are adjusted also destroys the
signature of MF. Therefore, it is almost impossible to observe the prominent
MF signature robust to the large-range external magnetic field and the gate
voltage in the huge mismatch of chemical potential.

Conclusion.—We define the dressed Majorana fermion (DMF) in any tunneling
strength and magnetic field for the hybrid nanowire. Under the open boundary,
the edge state of DMF is localized at both ends of the nanowire and
superconductor (SC). And we obtain the general topological phase of DMF, which
is determined by the magnetic field and the chemical potential in the nanowire
and SC. We clarify the validity of the low-energy phase diagram. Namely, it is
only an approximate result of the DMF phase in the weak tunneling scheme. We
also point out that in the mainstream hybrid InSb/As-Al system, the mismatch
in chemical potential between nanowire and SC makes MF challenging to observe
(Das  _et al._ , 2012; Zhang  _et al._ , 2018; Lutchyn  _et al._ , 2018;
Valentini  _et al._ , 2021).

Our theory provides a new method to define MF in the hybrid system
analytically. Namely, we treat the quasi-excitation in the nanowire and
superconductor equally instead of treating SC as an environment to provide the
proximity effect for the nanowire. Therefore, such a definition of DMF is
general and applicable to any parameters of any chosen nanowire and SC
materials, including the strong tunneling strength and magnetic field.

Note that although s-wave SC is described by a one-dimensional lattice model
in the above discussion, such a method of defining Majorana fermion can be
generalized to more hybrid systems, such as the two-dimensional multi-band
nanowire coupled to superconducting shell (Lutchyn  _et al._ , 2011; Stanescu
_et al._ , 2011; Reeg  _et al._ , 2017; Cole  _et al._ , 2016), and two-
dimensional topological insulator-superconductor system (Fu and Kane, 2008;
Wang  _et al._ , 2015; Yue  _et al._ , 2023).

The authors appreciate quite much for the helpful discussion with Sheng-Wen Li
in BIT, and Dong E. Liu, Zhan Cao and Gu Zhang in BAQIS. This study is
supported by the National Natural Science Foundation of China (NSFC) (Grant
No. 12088101) and NSAF (Grant No. U2230402).

## References

  * Kitaev (2001) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Sau  _et al._ (2010) J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Phys. Rev. Lett. 104, 040502 (2010).
  * Alicea (2012) J. Alicea, Reports on Progress in Physics 75, 076501 (2012).
  * Lutchyn  _et al._ (2018) R. M. Lutchyn, E. P. A. M. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, Nature Reviews Materials 3, 52 (2018).
  * Prada  _et al._ (2020) E. Prada, P. San-Jose, M. W. A. de Moor, A. Geresdi, E. J. H. Lee, J. Klinovaja, D. Loss, J. Nygård, R. Aguado, and L. P. Kouwenhoven, Nature Reviews Physics 2, 575 (2020).
  * Lutchyn  _et al._ (2010) R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg  _et al._ (2010) Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Mourik  _et al._ (2012) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
  * Stanescu  _et al._ (2011) T. D. Stanescu, R. M. Lutchyn, and S. Das Sarma, Phys. Rev. B 84, 144522 (2011).
  * Stanescu and Das Sarma (2017) T. D. Stanescu and S. Das Sarma, Phys. Rev. B 96, 014510 (2017).
  * (11) See supplemental material for details: The approximate conditions of low-energy effective theory and the phase diagram determined by it in Sec. I; The definition of the dressed Majorana fermion and the conditions under which it exists in Sec. II A; The topological phase of dressed Majorana fermion in Sec. II B ; Proof that the dressed Majorana fermion results in a zero-bias peak Sec. II C.
  * Qiao  _et al._ (2022) G.-J. Qiao, S.-W. Li, and C. P. Sun, Phys. Rev. B 106, 104517 (2022).
  * Das  _et al._ (2012) A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Nature Physics 8, 887 (2012).
  * Valentini  _et al._ (2021) M. Valentini, F. Peñaranda, A. Hofmann, M. Brauns, R. Hauschild, P. Krogstrup, P. San-Jose, E. Prada, R. Aguado, and G. Katsaros, Science 373, 82 (2021).
  * Zhang  _et al._ (2018) H. Zhang, C.-X. Liu, S. Gazibegovic, D. Xu, J. A. Logan, G. Wang, N. van Loo, J. D. S. Bommer, M. W. A. de Moor, D. Car, R. L. M. O. het Veld, P. J. van Veldhoven, S. Koelling, M. A. Verheijen, M. Pendharkar, D. J. Pennachio, B. Shojaei, J. S. Lee, C. J. Palmstrøm, E. P. A. M. Bakkers, S. D. Sarma, and L. P. Kouwenhoven, Nature 556, 74 (2018).
  * Vaitiekenas  _et al._ (2020) S. Vaitiekenas, G. W. Winkler, B. van Heck, T. Karzig, M.-T. Deng, K. Flensberg, L. I. Glazman, C. Nayak, P. Krogstrup, R. M. Lutchyn, and C. M. Marcus, Science 367, eaav3392 (2020).
  * Gazibegovic  _et al._ (2017) S. Gazibegovic, D. Car, H. Zhang, S. C. Balk, J. A. Logan, M. W. A. de Moor, M. C. Cassidy, R. Schmits, D. Xu, G. Wang, P. Krogstrup, R. L. M. O. het Veld, K. Zuo, Y. Vos, J. Shen, D. Bouman, B. Shojaei, D. Pennachio, J. S. Lee, P. J. van Veldhoven, S. Koelling, M. A. Verheijen, L. P. Kouwenhoven, C. J. Palmstrøm, and E. P. A. M. Bakkers, Nature 548, 434 (2017).
  * He  _et al._ (2017) Q. L. He, L. Pan, A. L. Stern, E. C. Burks, X. Che, G. Yin, J. Wang, B. Lian, Q. Zhou, E. S. Choi, K. Murata, X. Kou, Z. Chen, T. Nie, Q. Shao, Y. Fan, S.-C. Zhang, K. Liu, J. Xia, and K. L. Wang, Science 357, 294 (2017).
  * Stoudenmire  _et al._ (2011) E. M. Stoudenmire, J. Alicea, O. A. Starykh, and M. P. Fisher, Phys. Rev. B 84, 014503 (2011).
  * Chevallier and Klinovaja (2016) D. Chevallier and J. Klinovaja, Phys. Rev. B 94, 035417 (2016).
  * Cao  _et al._ (2022) Z. Cao, D. E. Liu, W.-X. He, X. Liu, K. He, and H. Zhang, Phys. Rev. B 105, 085424 (2022).
  * Jiang  _et al._ (2022) Y. Jiang, S. Yang, L. Li, W. Song, W. Miao, B. Tong, Z. Geng, Y. Gao, R. Li, F. Chen, Q. Zhang, F. Meng, L. Gu, K. Zhu, Y. Zang, R. Shang, Z. Cao, X. Feng, Q.-K. Xue, D. E. Liu, H. Zhang, and K. He, Phys. Rev. Mater. 6, 034205 (2022).
  * Reeg  _et al._ (2017) C. Reeg, D. Loss, and J. Klinovaja, Phys. Rev. B 96, 125426 (2017).
  * Clogston (1962) A. M. Clogston, Phys. Rev. Lett. 9, 266 (1962).
  * Maki and Tsuneto (1964) K. Maki and T. Tsuneto, Progress of Theoretical Physics 31, 945 (1964).
  * Reeg and Maslov (2017) C. Reeg and D. L. Maslov, Phys. Rev. B 95, 205439 (2017).
  * Vaitiekėnas  _et al._ (2018) S. Vaitiekėnas, M.-T. Deng, J. Nygård, P. Krogstrup, and C. M. Marcus, Phys. Rev. Lett. 121, 037703 (2018).
  * Li  _et al._ (2014) S.-W. Li, Z.-Z. Li, C. Y. Cai, and C. P. Sun, Phys. Rev. B 89, 134505 (2014).
  * Flensberg (2010) K. Flensberg, Phys. Rev. B 82, 180516 (2010).
  * Mikkelsen  _et al._ (2018) A. E. G. Mikkelsen, P. Kotetes, P. Krogstrup, and K. Flensberg, Phys. Rev. X 8, 031040 (2018).
  * Wang  _et al._ (2022) Z. Wang, H. Song, D. Pan, Z. Zhang, W. Miao, R. Li, Z. Cao, G. Zhang, L. Liu, L. Wen, R. Zhuo, D. E. Liu, K. He, R. Shang, J. Zhao, and H. Zhang, Phys. Rev. Lett. 129, 167702 (2022).
  * Liu  _et al._ (2017) C.-X. Liu, J. D. Sau, T. D. Stanescu, and S. Das Sarma, Phys. Rev. B 96, 075161 (2017).
  * Mohapatra  _et al._ (2019) S. Mohapatra, S. Mathimalar, S. Chaudhary, and K. Raman, Journal of Physics Communications 3, 045005 (2019).
  * Pan and Das Sarma (2020) H. Pan and S. Das Sarma, Phys. Rev. Res. 2, 013377 (2020).
  * Pan  _et al._ (2024) X.-H. Pan, L. Chen, D. E. Liu, F.-C. Zhang, and X. Liu, Phys. Rev. Lett. 132, 036602 (2024).
  * Lutchyn  _et al._ (2011) R. M. Lutchyn, T. D. Stanescu, and S. Das Sarma, Phys. Rev. Lett. 106, 127001 (2011).
  * Cole  _et al._ (2016) W. S. Cole, J. D. Sau, and S. Das Sarma, Phys. Rev. B 94, 140505 (2016).
  * Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Wang  _et al._ (2015) J. Wang, Q. Zhou, B. Lian, and S.-C. Zhang, Phys. Rev. B 92, 064520 (2015).
  * Yue  _et al._ (2023) X. Yue, G.-J. Qiao, and C. P. Sun, Phys. Rev. B 108, 195405 (2023).


# FMO Interaction Energy between 17$\beta$-Estradiol, 17$\alpha$-Estradiol and
Human Estrogen Receptor $\alpha$

Ricardo Ugarte∗ ∗Instituto de Ciencias Quimicas, Facultad de Ciencias,
Universidad Austral de Chile. Independencia 641, Valdivia, Chile (e-mail:
rugarte@uach.cl).

###### Abstract

The estrogen receptor is a nuclear hormone receptor activated by the natural
steroid hormone 17$\beta$-estradiol (E2). Fragment molecular orbital (FMO)
calculations were performed which allowed us to obtain the interaction energy
($E_{int}$) between E2, 17$\alpha$-estradiol (17$\alpha$-E2) and the human
estrogen receptor $\alpha$ ligand-binding domain. In aqueous media the
MP2/6-31G(d) $E_{int}$ was of -88.52 kcal/mol for E2 and -78.73 kcal/mol for
17$\alpha$-E2. Attractive dispersion interactions were observed between
ligands and all surrounding hydrophobic residues. Water molecules were found
at the binding site and strong attractive electrostatic interactions were
observed between the ligands and the Glu 353 and His 524 residues. The
essential dynamics revealed that E2 adapts to the binding site and its motion,
in a sense, synchronizes with the whole receptor; while 17$\alpha$-E2, with
its motion of greater amplitude compared to E2, disturbs the binding site.
Perhaps this feature of the normal substrate is a necessary condition for
biological function. Another important requirement relates to the number of
water molecules at the binding site. Therefore, negative values in $E_{int}$
is a necessary but not sufficient condition since, it is also necessary to
consider the conformers population that fulfill all the requirements that
ensure a biological response.

###### Index Terms:

Estrogen Receptor, Estradiol, FMO Calculations, Molecular Dynamics Simulation,
Clustering, Essential Dynamics, PCA Analysis

##  I Introduction

The estrogen receptor (ER) is a ligand-inducible transcription factor and a
member of the nuclear receptor superfamily that responds to endogenous
(endobiotics) and exogenous (xenobiotics) chemicals by regulating gene
expression [1]. ER include two subtypes, ER$\alpha$ and ER$\beta$. These share
structural patterns that are responsible for similar functional features. Both
subtypes possess a modular organization that is characteristic of the nuclear
receptors; five functional domains from the $NH_{2}-$ to $COOH$ terminus,
designated A/B, C (DNA-binding domain or DBD), D, E (ligand-binding domain or
LBD), and F [2, 3]. The full-length ER is not amenable to forming crystals
make x-ray crystallography impossible.

There are three main mammalian endobiotics called estrogens: estrone (E1),
17$\beta$-estradiol (E2), and estriol (E3), all derived from cholesterol. All
these estrogens act as ligands (L) for ER, but E2 (Figure 1) is the most
abundant and physiologically active form as it exhibits extremely high binding
affinity to ERs [4]. 17$\alpha$-estradiol (17$\alpha$-E2) is a stereoisomer of
E2 (C17 epimer of E2), and according to somewhat contradictory data it would
present an apparently lower estrogenic potency than E2 [5, 6, 7, 8]. The
xenobiotics that interfere with the normal ER signaling include
pharmaceuticals, industrial chemicals, pesticides, and phytoestrogens [2].
Estradiol binding induces a major structural reorganization of the LBD that
converts the inactive ER to the functionally active form by generating
surfaces for enhanced stability of the ER dimer and of the interacting
co-regulatory proteins [9].

Numerous crystal structures have been determined for the LBDs of both
subtypes, and these have given a detailed insight into the structure and
alterations during the ligand binding. The structure of ER LBD reveals a
conserved core of twelve $\alpha$-helices and a short two-stranded
antiparallel $\beta$-sheet arranged into a three-layered sandwich fold; this
arrangement generates a mostly hydrophobic cavity in the lower part of the
domain which can accommodate the ligand [1] (Figure 2).

A number of experimental and theoretical studies have been performed to
investigate the L-ER interaction [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34], and since 1997 about 100
crystal structures of ER LBD with several ligands have been solved and
deposited in the Research Collaboratory for Structural Bioinformatics (RCSB)
Protein Data Bank (PDB). On the basis of the above information, the mode of
binding between ERs and their ligands has been determined. The specific
recognition between ER and its ligand mainly depends on hydrogen bonds and
hydrophobic contacts [35, 36].

Most of the theoretical studies which use the structures deposited in RCSB PDB
are carried out using molecular dynamics (MD) simulations. These simulations
are based on empirical force fields that may not be accurate enough to predict
L-ER interaction energies. Accuracy requirements could be provided by ab
initio quantum mechanical calculations, but these can be very computationally
expensive and time consuming. The hybrid QM/MM (quantum mechanics/molecular
mechanics) is a method that combines the precision of quantum mechanics and
the speed of empirical force fields. In this approach, part of the system that
includes the chemically relevant region is treated quantum mechanically (QM)
while the remainder, often referred to as the environment, is treated at the
classical level using empirical or molecular mechanics (MM) force fields. This
multiscale approach reduces the computational cost significantly as compared
to a QM treatment of the entire system and makes simulations possible [37,
38].

Figure 1: Chemical structure of 17$\beta$-Estradiol (top) and
17$\alpha$-Estradiol (bottom). Figure 2: Model of HER$\alpha$ LBD (ribbon). E2
and water molecule (ball and stick) at the binding site [10]. The model based
on the RSCB PDB crystal structure (PDB code 1G50) includes 247 amino acid
residues.

An efficient alternative to either the full ab initio QM, MM or QM/MM, lies in
the fragment-based methods, which form an actively developed field of research
[39]. Fragment Molecular Orbital (FMO) method [40, 41] is one such method that
has been used for efficient and accurate QM calculations in very large
molecular systems [21, 42]. FMO involve fragmentation of the chemical system,
and from ab initio or density functional quantum-mechanical calculations of
each fragments (monomers) and their dimers (and trimers if greater accuracy is
required) one can construct the total properties. The distinctive feature of
FMO is the inclusion of the electrostatic field from the whole system in each
individual fragment calculation, and in using the systematic many-body
expansion. The FMO method is suited to various analyses, as it provides
information on fragments and their interactions that are naturally built into
the method.

In the present article, we report a study on an calculation of the interaction
energy between both estradiol isomers and LBD of human estrogen receptor alpha
in aqueous medium. The aim is to further explore the mechanism of interaction
underlying to the L-ER binding. In particular, verify if the calculation
method is capable of discriminating between two very similar compounds, which
bind to the same binding site. Briefly, the main steps of the calculation are
as follows: (i) Search for representative structures of the conformational
space around the crystallographic state of L-ER$\alpha$ LBD by means of MD
simulations and cluster analysis. (ii) Geometry optimization of the
representative structures using QM/MM approach. (iii) Single point FMO
calculation on the above optimized structures in order to obtain the inter-
fragment interaction energies.

##  II Methods

###  II-A Model Building

Crystal structure of the HER$\alpha$ LBD in complex with E2 (PDB code 1G50) at
2.9 Å resolution was retrieved from the RSCB PDB [43]. The model was built
from chain A of the homodimer. Missing hydrogen atoms were added with the LEaP
module of AmberTools 15 package [44]. AMBER FF14SB force field was selected
for the proteins and general AMBER force field (GAFF) parameters [45] were
employed for E2. In order to parameterize E2, electrostatic potential was
calculated by Gaussian 09 program at the HF/6-31G(d) level of theory [46].
Partial charges were fitted by RESP method of the Antechamber module of
AmberTools 15 [47]. Asp 313, Asp 321 and Glu 323 were modeled as neutral
species, Arg, Lys and remainder of Asp, Glu residues as charged species, all
tyrosines as neutral. Three of the 13 residues of histidine (His 476, His 501,
His 513), were protonated in order to preserve the electroneutrality of the
system. The N- and C-terminus residues were protonated and deprotonated,
respectively. The model was solvated with TIP3P water in a pre-equilibrated
box measuring 89 x 78 x 78 Å3. The E2-ER-W system contains 247 amino acid
residues (4009 atoms), E2 and 13394 water molecules (W). The total number of
atoms in the system is 44235.

Finally, the same previous protocol was used to build 17$\alpha$-E2-ER-W and
ER-W. The molecular structure of 17$\alpha$-E2 in the binding site was derived
from E2 by modifying the orientation of the corresponding OH group, while the
ER-W structure was generated by eliminating the ligand.

###  II-B Molecular Dynamics Simulations

In order to remove the steric strain introduced when adding hydrogen atoms,
holo (E2-ER-W, 17$\alpha$-E2-ER-W) and apo (ER-W) models were subjected to
three successive steps of minimization using the SANDER module of the
AmberTools 15. First, 1000 steps of steepest descent followed of 1500 steps of
conjugate gradient, allowing only H atoms and water to move while holding the
rest of the system fixed. Next, the same minimization algorithm is used as the
previous one, but maintaining the water fixed. Finally, the whole system was
minimized without any restraints for 2000 steps of steepest descent followed
by 1000 steps of conjugate gradient.

All simulations were carried out with the SANDER module of the AmberTools 15
with periodic boundary conditions, using Particle Mesh Ewald method [48] to
treat long-range electrostatics interactions with a non-bonded cutoff of 10 Å.
All bonds involving hydrogen atoms were restrained using the SHAKE [49]
algorithm. Temperature regulation was done using a Langevin thermostat with
collision frequency of 1 ps-1. The Berendsen barostat was used for constant
pressure simulation at 1 atm, with a relaxation time of 1 ps. The time step
was 1 or 2 fs. The hydrogen mass distribution (HMR) method was used for
accelerating MD simulations [50]. The final energy-minimized models was then
submitted to the following protocol:

Scheme 1:

NVT: 0 $\rightarrow$ 310 K $\Delta t=1$ fs 100 psNPT: 310 K $\Delta t=2$ fs
500 psNVT: 310 $\rightarrow$ 5 K $\Delta t=1$ fs 100 ps

From the restart file of the last simulation in Scheme 1, we conducted an
extensive set of molecular dynamics simulations to explore the conformational
space in the vicinity of the crystallographic structure. To circumvent the
limited conformational sampling ability of MD simulations at 310 K, we used
multiple-trajectory short-time simulations [51]. By combining the sampling
ability of the multiple trajectories, we expect to sample more conformational
space than single trajectory of the same length. The aforementioned restart
file was used as seed for 30 short-time simulations that obey the protocol
established in Scheme 2:

Scheme 2:

NVT: 5 $\rightarrow$ 150 K $\Delta t=2$ fs 60 psNPT: 150 $\rightarrow$ 310 K
$\Delta t=2$ fs 140 psNVE: 310 K $\Delta t=2$ fs 500 ps

The initial velocities (Scheme 2) were assigned randomly from a Maxwell-
Boltzmann distribution at 5 K. The trajectories start with the same structure
and differ only in the initial velocity assignment. At the end of the
equilibration, from $\sim$ 80 ps NPT ensemble, the average temperature of the
final 60 ps was 310 K, and the average density was 1.0 g/mL. All production
runs of 0.5 ns were performed in an NVE ensemble at 310 K.

###  II-C Cluster Analysis

Table I: Cluster Analysis of MD trajectories in the different models

| E2-ER-W | 17$\alpha$-E2-ER-W | ER-W
---|---|---|---
Cluster | 1 | 2 | 3 | 4 | 5 |  1 | 2 | 3 | 4 | 5 |  1 | 2 | 3 | 4 | 5
CP(a) | 3630 | 1749 | 1169 | 907 | 45 |  5156 | 185 | 543 | 77 | 1539 |  1800 | 3759 | 1652 | 191 | 98
CP % | 48.4 | 23.3 | 15.6 | 12.1 | 0.6 |  68.8 | 2.5 | 7.2 | 1.0 | 20.5 |  24.0 | 50.1 | 22.0 | 2.6 | 1.3
RS(b) |  $\beta$I |  $\beta$II |  $\beta$III |  $\beta$IV |  $\beta$V |  $\alpha$I |  $\alpha$II |  $\alpha$III |  $\alpha$IV |  $\alpha$V |  I | II | III | IV | V

  * •

(a)Cluster population. (b)Representative structure of the respective cluster.

Cluster analysis methods have been developed for analyzing simulation
trajectories of biomacromolecules and used for the analysis of their
conformational behavior in solution [52, 53]. These methods group together
similar conformers from molecular simulations. A clustering approach for each
model, based on the C$\alpha$-RMSD (root mean square deviation) was applied to
the snapshots of the MD simulations [54]. We selected the alpha carbon atoms
because they describe the backbone conformations. The C$\alpha$-RMSD was
calculated after rigid body alignment of C$\alpha$ atoms of each frame of the
trajectory with respect to C$\alpha$ atoms of the average structure of the
protein. Prior to clustering, the individual trajectories from the 30 short-
time simulations were combined into a single file (full trajectory), and the
water molecules were removed to speed up the calculations. In our analysis,
7500 snapshots were grouped into five clusters. Each cluster is described by a
centroid structure, which in itself is not physically significant as it is
effectively a mathematical construct based on the members of a cluster.
However, the actual structure closest to the centroid (rmsd) is significant
and representative of each cluster. Thus, five representative structures (RS)
of each cluster, and therefore of the conformers population, were obtained
(Table I).

###  II-D QM/MM Optimization

All L-ER-W models representative of the population were subjected to geometry
optimization using Gaussian 09 at the MM/AMBER level of theory. Then, to
facilitate the QM/MM calculations, the water molecules beyond 10 Å of the
protein surface were deleted using VMD program [55]. As a consequence, new
models (L-ER-w) with a water layer of 10 Å around of receptor surface were
generated. ONIOM, [37] a hybrid QM/MM method implemented in Gaussian 09, was
used for the geometry optimization of L-ER-w models. In the present study we
used a two-layer ONIOM (B3LYP/6-31G(d):AMBER) scheme: L(B3LYP/6-31G(d));
ER-w(AMBER).

###  II-E FMO Calculations

Fragment-based approaches refer to the chemical idea of parts of the system
retaining their identity to a large extent (e.g., functional groups and
residues). FMO method not only reduces the computational cost, but it also
provides a wealth of information on the properties of fragments and their
interactions. The calculation of the ligand-receptor interaction energy is
based on obtaining and sum the interaction energies between all pairs of
fragments ligand-amino acid residue and if applicable, ligand-water trapped on
binding site.

Table II: Number of FMO fragments in the L-ER-w models(a)

Fragments | E2-ER-w | 17$\alpha$-E2-ER-w
---|---|---
17$\beta$-estradiol | 1 |
17$\alpha$-estradiol |  | 1
Amino acid residue | 246 | 246
Water | 4179 | 4179
Total | 4426 | 4426

  * •

(a)Models with a water layer of 10 Å around of receptor surface.

The AFO (adaptive frozen orbitals) scheme was used throughout for
fragmentation across peptide bonds, with the default settings for bond
definitions. The fragmentation of the model was as follows: The first two
amino acid residues and each remaining amino acid residue of apo-ER, L, and
the water molecule were treated as a single fragment. Table II shows the
number of fragments in the different models. In order to reduce the
computational cost, multilayer two-body FMO calculation were performed. The
molecular system was divided into two layers treated at different levels of
theory: FMO2-RHF/STO-3G:MP2/6-31G(d). The layer 1 (aqueous environment)
described by RHF/STO-3G and layer 2 (L-ER) described by MP2/6-31G(d) [56]. The
water molecules of the binding site were included in the layer 2.

To contrast the interaction energies, models in vacuo were built from the
original models in which all water molecules were removed, except those found
on the binding site; thus, these water molecules are incorporated into the
receptor. The E2-ER and 17$\alpha$-E2-ER were subjected to FMO2-MP2/6-31G(d)
calculations.

Finally, pair interaction energies were computed in the following models:
E2-ER-w, 17$\alpha$-E2-ER-w, E2-ER and 17$\alpha$-E2-ER.

##  III Results and Discussion

###  III-A Analysis of Molecular Dynamics Trajectories

The RMSD with respect to average structure or RMSF were calculated (Å) for all
three models: apo-ER: 0.76 $\pm$ 0.12; E2-ER: 0.75 $\pm$ 0.10;
17$\alpha$-E2-ER: 0.80 $\pm$ 0.12, which indicates that the simulations
represent fluctuations around a stable average. The RMSD value of the holo-ER
average structures with respect to the apo-ER average structure taken as
reference are 0.98 Å and 0.72 Å for E2-ER and 17$\alpha$-E2-ER, respectively
(Figure 3). Apparently 17$\beta$-estradiol affects the apo-ER backbone more
than 17$\alpha$-estradiol. Detailed observation of the structures shows a
slight increase of random coil in holo-ER with respect to apo-ER (the increase
is higher in E2-ER). The above, is in detriment of the content of alpha helix.
The following motifs are not influenced by the ligands binding: H1/H2, Coil,
H3, H6 and H10/H11. The loop H3-H4 is only affected by 17$\alpha$-E2 binding;
the loop H8-H9 is significantly affected by E2, and only slightly affected by
17$\alpha$-E2. The N- and C-terminus are also slightly modified by ligands
binding, although to a larger extent by E2. The $\beta$-hairpins is modified
in the same manner by the ligands binding. With regard to the ligands, both
are placed in the same position on the binding site (Figure 4).

Figure 3: Stereo view of the average structures over the trajectories: apo-ER
(tan); E2-ER (sky blue) and 17$\alpha$-E2-ER (pale plum). Estrogen Receptor:
ribbon; Ligands: ball and stick. Structure superposition based on the
C$\alpha$. Taking as a reference the apo-ER structure, the root mean square
distance (RMSD) value for E2-ER and 17$\alpha$-E2-ER is 0.98 Å and 0.72 Å,
respectively. Figure 4: 17$\beta$-Estradiol and 17$\alpha$-Estradiol on the
binding site. Detail from image taken from the Figure 3.

In order to compare the dynamic behavior of the three models, the essential
dynamics method was used [57, 58]. For this purpose, the following procedure
was carried out: (i) Combination of the full trajectories (water stripped)
into a single trajectory (22500 frames) and removal of ligand molecules. (ii)
Principal component analysis (PCA) [59] on the combined MD trajectory to
calculate a set of eigenvectors, and then projected the backbone C$\alpha$
fluctuations onto these eigenvectors or principal components. (iii) Clustering
of the combined MD trajectory and PCA data sorted by cluster [54].

Figure 5: Clustering results of MD simulations projected on 2D plane formed by
first two PCs. Clustering was performed on the combined MD trajectory data and
two clusters were requested in the computation. Cluster separation for the
backbone dynamics of the entire estrogen receptor (top left) and its binding
site (bottom left): cluster 1ra (purple); cluster 2r (orange); cluster 1bsb
(magenta); cluster 2bs (cyan). Distribution of motion into the clusters (top
and bottom right): ER (yellow); E2-ER (red); 17$\alpha$-E2-ER (blue); E2-ER
(diamond dark green) and 17$\alpha$-E2-ER (diamond dark pink) representative
structures. ar: receptor; bbs: binding site. Figure 6: Superposition of
representative structures from cluster 1r (tan) and 2r (sky blue). The root
mean square distance (RMSD) value for this superposition is 1.045 Å. Figure 7:
Probability distributions for the backbone C$\alpha$ fluctuations in combined
MD trajectory along several principal components. The entire estrogen receptor
(top) and its binding site (bottom): apo-ER (yellow); E2-ER (red);
17$\alpha$-E2-ER (blue).

While PCA reveals the essential dynamics or main motions (large-scale motion)
contained in an MD trajectory it does not partition the frames into distinct
conformational categories. This can be achieved by clustering the PC data. A
previous study confirms the choice of two clusters. If three clusters are
requested, the third represents a population percentage of approximately 0.5%.
The trajectory of the atomic displacements or fluctuations described by
principal components (PCs) 1 and 2 are shown in Figure 5. Clustering for the
C$\alpha$ fluctuations of the entire estrogen receptor is displayed in figure
5 (top left). This trajectory contains two conformational categories,
practically separated along the first PC. Distribution of motion of apo-ER and
holo-ER into the clusters are shown (top right). Interestingly, the apo-ER
form displays dynamics characteristics of both clusters. The representative
structure of cluster 1r belongs to the E2-ER population, while that of cluster
2r belongs to the 17$\alpha$-E2-ER population (Figure 6). The motions along
these PCs are qualitatively similars since the point cloud represented in the
plane has analogous sizes and shapes for apo-ER and holo-ER forms. This seems
to indicate that there is no significant change in global dynamics of estrogen
receptor upon uptake of ligand. However, the simulations visit different areas
of the essential 2D space, indicating that there may be a small difference in
the conformations associated with each model. Remarkably, a subset of
trajectories of apo-ER and 17$\alpha$-E2-ER overlap; apparently, the overall
differences between these forms are not as pronounced.

Table III: FMO Interaction Energy

In Aqueous FMO2-RHF/STO-3G:MP2/6-31G(d)
---
| E2-ER-w | 17$\alpha$-E2-ER-w
RS |  $\beta$I |  $\beta$II |  $\beta$III |  $\beta$IV |  $\beta$V |  $\alpha$I |  $\alpha$II |  $\alpha$III |  $\alpha$IV |  $\alpha$V
$E_{L-Aa}$ | -88.47 | -82.82 | -82.91 | -82.43 | -79.37 | -77.54 | -81.36 | -75.18 | -67.44 | -81.73
$E_{L-w}$ | -3.21 | -0.70 | -6.90 | -1.50 | -7.70 |  | -8.69 | -6.62 | -11.43 | 1.45
$E_{int-RS}^{(a)}$ | -91.68 | -83.52 | -89.81 | -83.93 | -87.07 | -77.54 | -90.05 | -81.80 | -78.87 | -80.28
$E_{int}^{(b)}$ | -88.52 | -78.73
In Vacuo FMO2-MP2/6-31G(d)
| E2-ER | 17$\alpha$-E2-ER
RS |  $\beta$I |  $\beta$II |  $\beta$III |  $\beta$IV |  $\beta$V |  $\alpha$I |  $\alpha$II |  $\alpha$III |  $\alpha$IV |  $\alpha$V
$E_{L-Aa}$ | -91.80 | -86.34 | -87.22 | -85.47 | -80.78 | -80.93 | -82.84 | -79.94 | -69.16 | -84.40
$E_{L-w}$ | -2.829 | -0.755 | -6.256 | -1.140 | -7.925 |  | -8.784 | -6.110 | -11.78 | 1.37
$E_{int-RS}^{(a)}$ | -94.63 | -87.10 | -93.48 | -86.61 | -88.70 | -80.93 | -91.62 | -86.05 | -80.94 | -83.03
$E_{int}^{(b)}$ | -91.69 | -82.00

  * •

$L$: Ligand; $Aa$: Amino acids; $w$: waters. (a) $E_{int-
RS}=(E_{L-Aa}+E_{L-w})$. (b)Represents the ligand-receptor interaction energy
and correspond to the average weighted energy based on the cluster population.
All energies at MP2/6-31G(d) level of theory in kcal/mol.

Table IV: FMO Interaction Energy

In Aqueous FMO2-RHF/STO-3G:MP2/6-31G(d)
---
| E2-ER-w | 17$\alpha$-E2-ER-w
RS |  $\beta$I |  $\beta$II |  $\beta$III |  $\beta$IV |  $\beta$V |  $\alpha$I |  $\alpha$II |  $\alpha$III |  $\alpha$IV |  $\alpha$V
$E_{L-Glu\ 353}$ | -24.51 | -24.50 | -24.21 | -23.16 | -22.92 | -22.55 | -27.37 | -21.00 | -2.24 | -24.92
$E_{L-Arg\ 394}$ | 3.24 | 4.40 | 2.45 | 2.56 | 2.55 | 1.01 | 1.04 | 1.76 | -9.07 | -0.06
$E_{L-His\ 524}$ | -16.82 | -16.42 | -16.45 | -15.60 | -16.61 | -15.86 | -13.09 | -11.96 | -14.01 | -13.41
In Vacuo FMO2-MP2/6-31G(d)
| E2-ER | 17$\alpha$-E2-ER
RS |  $\beta$I |  $\beta$II |  $\beta$III |  $\beta$IV |  $\beta$V |  $\alpha$I |  $\alpha$II |  $\alpha$III |  $\alpha$IV |  $\alpha$V
$E_{L-Glu\ 353}$ | -26.73 | -27.16 | -27.50 | -25.75 | -23.49 | -24.56 | -28.16 | -25.48 | -2.64 | -26.60
$E_{L-Arg\ 394}$ | 3.51 | 4.61 | 2.90 | 2.87 | 2.44 | 0.73 | 0.98 | 2.39 | -9.87 | -0.28
$E_{L-His\ 524}$ | -17.59 | -16.84 | -17.16 | -15.88 | -16.86 | -16.35 | -13.24 | -12.21 | -14.04 | -13.94

  * •

$L$: Ligand. All energies at MP2/6-31G(d) level of theory in kcal/mol.

In Figure 5 (bottom), cluster separation for the C$\alpha$ fluctuations of the
binding site is displayed. Cluster 1bs is mainly occupied by E2-ER, while
cluster 2bs by apo-ER; 17$\alpha$-E2-ER is distributed evenly between both
clusters. Obviously, the essential 2D space is smaller because the motions are
more restricted than those of entire receptor. By observing the dynamic
behavior and the small conformational differences in the models, it is
possible to argue that the most active ligand (e.g. 17$\beta$-Estradiol)
causes significant changes in the apo-receptor that eventually are translated
into biochemical responses. The above argument would explain the fact that
17$\alpha$-E2-ER exhibits, in the binding site, features of both ”antipode
structures” represented by E2-ER and apo-ER, and is distributed between both
clusters Figure 5 (bottom right).

Probability distributions for the backbone C$\alpha$ fluctuations in combined
MD trajectory along several principal components are shown in Figure 7. In
order to provide insights about the three PCA modes of motion, movies to
visualize pseudo-trajectories of the models (download source) were created
with the ambertools cpptraj module. The comparison of the top right panel in
figure 7 with the third PCA mode movies, taking into account the overlap of
the distributions, allows us to ignore the particularities of the motion
displayed in the movies and focus on its global features. Consequently,
comparing the PCA modes of motion and the respective probability distributions
we could conclude the following: (i) A greater plasticity in the movement of
apo-ER with respect to both holo-ER. The backbone rigidity in movement is
induced by the ligand (second PCA mode movies vs. top middle panel). (ii) A
wiggling motion through the principal axis of the estrogen receptor that
occurs at apo-ER and 17$\alpha$-E2-ER. This movement has a slightly smaller
amplitude in apo-ER (first PCA mode movies vs. top left panel). An analogous
analysis at the binding site, including the respective PCA mode movies, showed
two types of movement of the ligand; a jiggling motion and a seesaw motion
around a principal axis lying in the plane of the molecule. The amplitude of
motion was larger in 17$\alpha$-E2 than in E2, which results in a significant
disturbance of the binding site. In fact, the interaction with His 524 is
continuously hindered by the strong movement of 17$\alpha$-E2. The normal
substrate E2, fits perfectly to the binding site and its movement, to some
extent, harmonizes with the whole receptor. In general, when comparing the
models we can affirm that apo form show a greater plasticity in its movement,
and that E2-ER has less conformational freedom.

Figure 8: Stereo view of estrogen receptor binding site: 17$\beta$-estradiol
(ball and stick) is shown surrounded by a hydrophobic pocket (blue wire),
water molecule (ball and stick) and important amino acid residues (stick): His
524, Glu 353 and Arg 394. H-bond (cyan). The image corresponds to $\beta$I
representative structure. Figure 9: Stereo view of estrogen receptor binding
site: 17$\alpha$-estradiol (ball and stick) is shown surrounded by a
hydrophobic pocket (blue wire) and important amino acid residues (stick): His
524, Glu 353 and Arg 394. H-bond (cyan). The image corresponds to $\alpha$I
representative structure. Figure 10: Stereo view of estrogen receptor binding
site: 17$\alpha$-estradiol (ball and stick) is shown surrounded by a
hydrophobic pocket (blue wire), water molecule (ball and stick) and important
amino acid residues (stick): His 524, Glu 353 and Arg 394. H-bond (cyan). The
image corresponds to $\alpha$III representative structure.

###  III-B FMO Interaction Energies

From MD simulations, 7500 snapshots of each model were grouped into five
clusters. Thus, five representative structures (RS) of each cluster were
obtained. The interaction energies at MP2/6-31G(d) are shown in Table III. The
L-ER interaction energy ($E_{int}$) values are negative which, at least at
this level of calculation, would indicate a certain stability of the systems.

In a FMO study of HER$\alpha$ LBD in complex with 17$\beta$-estradiol (PDB
code 1ERE), the model included 241 amino acid residues, one water molecule
which directly mediates E2-ER binding (where the hydrogen bonded water
molecule was incorporated in the receptor) and E2. The binding energy was
estimated from: $\Delta E_{Binding}=E_{E2-ER}-(E_{ER}+E_{E2})$. The total
energies (E) were considered in the calculation and the geometries of ER and
E2 were fixed in those found in E2-ER model. The reported binding energy was
-37.65 kcal/mol at FMO2-RHF/STO-3G level of theory [21]. In another FMO study
with the same model system, FMO2 interaction energy between E2 and ER was
calculated using HF and MP2 methods with several basis sets [60]. The
calculated interaction energy was -40.26 kcal/mol at FMO2-RHF/6-31G(d) and
-123.73 kcal/mol at FMO2-MP2/6-31G(d) level of theory. The large interaction
energy difference between the RHF and the MP2 methods is due to dispersion
interaction, which can only be described by electron correlation methods. In
general, charged and polarized amino acid residues interact either strongly or
weakly with the ligand, while hydrophobic residues contribute to weak
interactions. The sum of these weak dispersion interactions makes the
difference between both methods. In the present study, interaction energy is
reported at the MP2/6-31G(d) level. The results of the calculation of the
lower level are excluded.

Table V: Number of water molecules in the binding site E2-ER-w |  $\beta$I |  $\beta$II |  $\beta$III |  $\beta$IV |  $\beta$V
---|---|---|---|---|---
water | 1 | 1 | 1 | 1 | 3
17$\alpha$-E2-ER-w |  $\alpha$I |  $\alpha$II |  $\alpha$III |  $\alpha$IV |  $\alpha$V
water | 0 | 2 | 1 | 2 | 1
ER-w | I | II | III | IV | V
water | 7 | 8 | 6 | 3 | 2

In holo models we observe that Glu 353 and His 524 present a strong
electrostatic interaction with the ligand (Table IV), and it is known that
together with Arg 394, they form a hydrogen bond network with E2. A water
molecule is similarly responsible for yet another stabilizing hydrogen bond
between E2 and the apo-ER [22, 60]. Many binding site hydrophobic residues are
stabilized (attractive interaction) through dispersion interactions with the
ligand; the MP2 electron correlation is essential to characterize these
interactions, whose function is possibly to accommodate the substrate.
Relevant hydrophobic residues at the binding site are: Leu 346, Leu 387, Leu
391, Phe 404 and Leu 525 [61]. Therefore, the interactions between the ligand
with these functionally important molecules, together with the hydrogen bonds
network play a key role in the L-ER binding (Figures 8-10).

During the simulation water molecules were trapped at the binding site (Table
V), which consists of all residues that have at least one atom within 4 Å from
any ligand atom in L-ER-w models. This generally gives a good representation
of the important residues in the binding pocket of a protein. The number of
water molecules varies according to the representative structure and the
L-water interactions can be either attractive or repulsive. Focusing on the
more populated clusters, which in general terms determine the biological
response, the $\alpha$I representative structure lacks a water molecule in the
binding site and possibly does not fulfill a basic condition for an adequate
attachment of the ligand, and therefore for a suitable biological function.
While $\alpha$III contains a water molecule, it represents only 7.2% of the
population of 17$\alpha$-E2-ER. The binding site of $\alpha$V also contains a
water molecule, but the phenolic A ring hydroxyl group points in another
direction, and probably is not a biologically viable structure. This aspect is
interesting since it is possible to infer that not all the structures that
present affinity could present efficacy, and therefore, would not be relevant
to biological function. That is, only a subset of the conformers population
should be biologically functional.

Evidently, the problem arises from the 17$\alpha$-hydroxyl group-His 524
interaction and propagates to the rest of the system. By pointing the
$C_{17}$-OH group (D ring) in a direction contrary to that of the normal
substrate, it disturbs the anchorage of the molecule in the binding site. This
is supported by previous discussion of the PCA modes of motion and justified,
to some extent, by the fact that the 17$\alpha$-E2-His 524 interaction energy
is on average $\approx$ 1.5 kcal/mol lower than E2-His 524 interaction energy.

Even though both ligands display affinity for the receptor (E2 in greater
measure than 17$\alpha$-E2), the population of probably biologically viable
conformations is lower in 17$\alpha$-E2 than in E2, and this should be
considered when speculating theoretically about the possible action of a
ligand on a receptor.

##  IV Conclusions

The molecular interactions between 17$\beta$-estradiol, 17$\alpha$-estradiol
and human estrogen receptor $\alpha$ were calculated, and from these the L-ER
interaction energy ($E_{int}$) was obtained. The FMO2-MP2/6-31G(d) $E_{int}$
was of -88.52 kcal/mol for E2-ER-w and -78.73 kcal/mol for 17$\alpha$-E2-ER-w
models. In vacuo, $E_{int}$ was of -91.69 kcal/mol and -82.00 kcal/mol for
E2-ER and 17$\alpha$-E2-ER models, respectively. The results obtained for
17$\alpha$-estradiol clearly suggests that it tends to form a complex with the
receptor.

In general, attractive dispersion interactions were observed between ligands
and hydrophobic residues; these interactions play an important role in
stabilizing E2 and 17$\alpha$-E2 at the binding site. Water molecules were
found at the binding site of all representative structures, except $\alpha$I.
Strong attractive electrostatic interactions were observed between the ligands
and the following charged/polarized residues: Glu 353, His 524. These residues
tend to be located at the ends of the ligands, close to the OH groups of the A
and D rings.

The C17 epimer of E2 disturbs the binding site due to its strong motion, while
E2 fits perfectly to the binding site and its motion of lower amplitude, in a
sense, synchronizes with the whole receptor. Perhaps this feature of the
normal substrate is a necessary condition for biological function. Another
important requirement relates to the number of water molecules at the binding
site. Therefore, negative values in $E_{int}$ is a necessary but not
sufficient condition since, it is also necessary to consider the conformers
population that fulfill all the requirements that ensure a biological
response.

## Acknowledgements

I wish to express my gratitude to arXiv repository and to each and every
person who has contributed to the development and maintenance of free and open
source software.

## References

  * [1] S. Lorenzetti, L. Narciso, Nuclear Receptors: Connecting Human Health to the Environment, in: P. Cozzini, G.E Kellogg (Eds.), Computational Approaches to Nuclear Receptors, RSC Publishing, Cambridge, 2012, pp. 1-22.
  * [2] I. Tsakovska, I Pajeva, P. Alov, A. Worth, Recent Advances in the Molecular Modeling of Estrogen Receptor-Mediated Toxicity,  Advances in Protein Chemistry and Structural Biology 85 (2011) 217-251. https://doi.org/10.1016/B978-0-12-386485-7.00006-5.
  * [3] S.C. Hewitt, K.S. Korach, Estrogen Receptors: New Directions in the New Millennium,  Endocrine Reviews 39 (2018) 664-675. https://doi.org/10.1210/er.2018-00087.
  * [4] S. Srinivasan, Z. Nawaz, Molecular Biology of Estrogen Receptor Action, in: D.W. Pfaff, A.P. Arnold, … R.T. Rubin (Eds.), Hormones, Brain and Behavior (Second Edition), Academic Press, 2009, pp. 1187-1218.
  * [5] B.A. Littlefieldf, E. Gurpide, L. Markiewicz, B. McKinley, R.B. Hochberg, A Simple and Sensitive Microtiter Plate Estrogen Bioassay Based on Stimulation of Alkaline Phosphatase in Ishikawa Cells: Estrogenic Action of $\Delta$5 Adrenal Steroids,  Endocrinology 127 (1990) 2757-2762. https://doi.org/10.1210/endo-127-6-2757.
  * [6] G.M. Anstead, K.E. Carlson, J.A. Katzenellenbogen, The estradiol pharmacophore: Ligand structure-estrogen receptor binding affinity relationships and a model for the receptor binding site,  Steroids 62 (1997) 268-303. https://doi.org/10.1016/S0039-128X(96)00242-5.
  * [7] B.T. Zhu, G-Z. Han, J-Y. Shim, Y. Wen, X-R. Jiang, Quantitative Structure-Activity Relationship of Various Endogenous Estrogen Metabolites for Human Estrogen Receptor $\alpha$ and $\beta$ Subtypes: Insights into the Structural Determinants Favoring a Differential Subtype Binding,  Endocrinology 147 (2006) 4132-4150. https://doi.org/10.1210/en.2006-0113.
  * [8] C.D. Toran-Allerand, A.A. Tinnikov, R.J. Singh, I.S. Nethrapalli, 17 $\alpha$-Estradiol: A Brain-Active Estrogen?,  Endocrinology 146 (2005) 3843-3850. https://doi.org/10.1210/en.2004-1616.
  * [9] P. Yasar, G. Ayaz, S.D User, G. Gupur, M. Muyan, Molecular mechanism of estrogen–estrogen receptor signaling,  Reprod. Med. Biol. 16 (2017) 4-20. https://doi.org/10.1002/rmb2.12006.
  * [10] E.F. Pettersen, T.D. Goddard, C.C. Huang, G.S. Couch, D.M. Greenblatt, E.C. Meng, T.E. Ferrin, UCSF Chimera–a visualization system for exploratory research and analysis,  J. Comput. Chem. 25 (2004) 1605-1612. https://doi.org/10.1002/jcc.20084.
  * [11] S. Sasson, A.C. Notides, Estriol and Estrone Interaction with the Estrogen Receptor. II. Estriol and estrone-induced inhibition of the cooperative binding of [3H]estradiol to the estrogen receptor,  J. Biol. Chem. 258 (1983) 8118-8122.
  * [12] B.S. Katzenellenbogen, Biology and receptor interactions of estriol and estriol derivatives in vitro and in vivo,  J. Steroid. Biochem. 20 (1984) 1033-1037.
  * [13] G.G.J.M. Kuiper, J.G. Lemmen, B. Carlsson, J.C. Corton, S.H. Safe, P.T. van der Saag, B. van der Burg, J-A. Gustafsson, Interaction of Estrogenic Chemicals and Phytoestrogens with Estrogen Receptor $\beta$,  Endocrinology 139 (1998) 4252-4263. https://doi.org/10.1210/endo.139.10.6216.
  * [14] K.W. Gaido, S.C. Maness, D.P. McDonnell, S.S. Dehal, D. Kupfer, S. Safe, Interaction of Methoxychlor and Related Compounds with Estrogen Receptor alpha and beta, and Androgen Receptor: Structure-Activity Studies,  Mol. Pharmacol. 58 (2000) 852-858. https://doi.org/10.1124/mol.58.4.852.
  * [15] A. Tamrazi, K.E. Carlson, J.A. Katzenellenbogen, Molecular Sensors of Estrogen Receptor Conformations and Dynamics,  Mol. Endocrinol. 17 (2003) 2593-602. https://doi.org/10.1210/me.2003-0239.
  * [16] A. Matsushima, X. Liu, H. Okada, M. Shimohigashi, Y. Shimohigashi, Bisphenol AF Is a Full Agonist for the Estrogen Receptor ER$\alpha$ but a Highly Specific Antagonist for ER$\beta$,  Environmental Health Perspectives 118 (2010) 1267-1272. https://dx.doi.org/10.1289
  * [17] C.N. Harvey, J.C. Chen, C.A. Bagnell, M. Uzumcu, Methoxychlor and Its Metabolite HPTE Inhibit cAMP Production and Expression of Estrogen Receptors $\alpha$ and $\beta$ in the Rat Granulosa Cell In Vitro,  Reprod Toxicol. 51 (2015) 72-78. https://dx.doi.org/10.1016
  * [18] H. Cao, F. Wang, Y. Liang, H. Wang, A. Zhang, M. Song, Experimental and computational insights on the recognition mechanism between the estrogen receptor $\alpha$ with bisphenol compounds,  Arch. Toxicol. 91 (2017) 3897-3912. https://doi.org/10.1007/s00204-017-2011-0.
  * [19] R. Mesnage, A. Phedonos, M. Biserni, M. Arno, S. Balu, J.C. Corton, R. Ugarte, M.N. Antoniou, Evaluation of estrogen receptor alpha activation by glyphosate-based herbicide constituents,  Food and Chem. Toxicol. 108 (2017) 30-42. https://doi.org/10.1016/j.fct.2017.07.025.
  * [20] M.M.H. van Lipzig, A.M. ter Laak, A. Jongejan, N.P.E. Vermeulen, M. Wamelink, D. Geerke, J.H.N. Meerman, Prediction of Ligand Binding Affinity and Orientation of Xenoestrogens to the Estrogen Receptor by Molecular Dynamics Simulations and the Linear Interaction Energy Method,  J. Med. Chem. 47 (2004) 1018-1030. https://doi.org/10.1021/jm0309607.
  * [21] K. Fukuzawa, K. Kitaura, M. Uebayasi, K. Nakata, T. Kaminuma, T. Nakano, Ab initio Quantum Mechanical Study of the Binding Energies of Human Estrogen Receptor $\alpha$ with Its Ligands: An Application of Fragment Molecular Orbital Method,  J. Comput. Chem. 26 (2005) 1-10. https://doi.org/10.1002/jcc.20130.
  * [22] L. Celik, J.D.D. Lund, B. Schiott, Conformational Dynamics of the Estrogen Receptor $\alpha$: Molecular Dynamics Simulations of the Influence of Binding Site Structure on Protein Dynamics,  Biochemistry 46 (2007) 1743-1758. https://doi.org/10.1021/bi061656t.
  * [23] M.T. Sonoda, L. Martinez, P. Webb, M.S. Skaf, I. Polikarpov, Ligand Dissociation from Estrogen Receptor Is Mediated by Receptor Dimerization: Evidence from Molecular Dynamics Simulations,  Mol. Endocrinol. 22 (2008) 1565–1578. https://dx.doi.org/10.1210
  * [24] T.D. McGee, J. Edwards, A.E. Roitberg, Preliminary Molecular Dynamic Simulations of the Estrogen Receptor Alpha Ligand Binding Domain from Antagonist to Apo, Int. J.  Environ. Res. Public. Health 5 (2008) 111-114.
  * [25] F. Spyrakis, P. Cozzini, How Computational Methods Try to Disclose the Estrogen Receptor Secrecy-Modeling the Flexibility,  Curr. Med. Chem. 16 (2009) 2987-3027. https://doi.org/10.2174/092986709788803123.
  * [26] S. Zhuang, J. Zhang, Y. Wen, C. Zhang, W. Liu, Distinct mechanisms of endocrine disruption of DDT-related pesticides toward estrogen receptor $\alpha$ and estrogen-related receptor $\gamma$,  Environ. Toxicol. Chem. 31 (2012) 2597-605. https://doi.org/10.1002/etc.1986.
  * [27] L. Gao, Y. Tu, H. Ågren, L.A. Eriksson, Characterization of Agonist Binding to His524 in the Estrogen Receptor $\alpha$ Ligand Binding Domain,  J. Phys. Chem. B 116 (2012) 4823-4830. https://doi.org/10.1021/jp300895g.
  * [28] D. Jereva, F. Fratev, I. Tsakovska, P. Alov, T. Pencheva, I. Pajeva, Molecular dynamics simulation of the human estrogen receptor alpha: contribution to the pharmacophore of the agonists,  Mathematics and Computers in Simulation 133 (2017) 124-134. https://doi.org/10.1016/j.matcom.2015.07.003.
  * [29] L. Li, Q. Wang, Y. Zhang, Y. Niu, X. Yao, H. Liu, 2015. The Molecular Mechanism of Bisphenol A (BPA) as an Endocrine Disruptor by Interacting with Nuclear Receptors: Insights from Molecular Dynamics (MD) Simulations.  PLoS One. 10, e0120330. https://dx.doi.org/10.1371
  * [30] A. Zafar, S. Ahmadb, I. Naseem, 2015. Insight into the structural stability of Coumestrol with Human Estrogen Receptor $\alpha$ and $\beta$ subtypes: A combined approach involving docking and molecular dynamics simulation studies.  Wiley RSC Adv. 5, 81295. https://doi.org/10.1039/C5RA14745J.
  * [31] A. Shtaiwi, R. Adnan, M. Khairuddean, M. Al-Qattan, 2018. Molecular dynamics simulation of human estrogen receptor free and bound to morpholine ether benzophenone inhibitor.  Theor. Chem. Acc. 137, 101. https://doi.org/10.1007/s00214-018-2277-1.
  * [32] M. Pavlin, A. Spinello, M. Pennati, N. Zaffaroni, S. Gobbi, A. Bisi, G. Colombo, A. Magistrato, 2018. A Computational Assay of Estrogen Receptor $\alpha$ Antagonists Reveals the Key Common Structural Traits of Drugs Effectively Fighting Refractory Breast Cancers.  Sci. Rep. 8, 649. https://doi.org/10.1038/s41598-017-17364-4.
  * [33] E.E.M. Eid, F. Azam, M. Hassan, I. M. Taban, M.A. Halim, Zerumbone binding to estrogen receptors: an in-silico investigation,  Journal of Receptors and Signal Transduction 38 (2018) 342-351. https://doi.org/10.1080/10799893.2018.1531886.
  * [34] H. Cao, L. Wang, M. Cao, T. Ye, Y. Sun, Computational insights on agonist and antagonist mechanisms of estrogen receptor $\alpha$ induced by bisphenol A analogues,  Environmental Pollution 248 (2019) 536-545. https://doi.org/10.1016/j.envpol.2019.02.058.
  * [35] K. Kato, K. Fujii, T. Nakayoshi, Y. Watanabe, S. Fukuyoshi, K. Ohta, Y. Endo, N. Yamaotsu, S. Hirono, E. Kurimoto, A. Oda, 2018. Structural differences between the ligand-binding pockets of estrogen receptors alpha and beta.  J. Phys.: Conf. Ser. 1136, 012021. https://doi.org/10.1088/1742-6596/1136/1/012021.
  * [36] S. Lee, M.G. Barron, 2017. Structure-Based Understanding of Binding Affinity and Mode of Estrogen Receptor $\alpha$ Agonists and Antagonists.  PLoS One. 12, e0169607. https://dx.doi.org/10.1371
  * [37] T. Vreven, K. Morokuma, Hybrid Methods: ONIOM(QM:MM) and QM/MM,  Annual Reports in Computational Chemistry 2 (2006) 35-51. https://doi.org/10.1016/S1574-1400(06)02003-2.
  * [38] A.W. Götz, M.A. Clark, R.C. Walker, An Extensible Interface for QM/MM Molecular Dynamics Simulations with AMBER,  Journal of Computational Chemistry 35 (2014) 95-108. https://doi.org/10.1002/jcc.23444.
  * [39] M.S. Gordon, D.G. Fedorov, S.R. Pruitt, L.V. Slipchenko, Fragmentation Methods: A Route to Accurate Calculations on Large Systems,  Chem. Rev. 112 (2012) 1632-672. https://doi.org/10.1021/cr200093j.
  * [40] D. Fedorov, K. Kitaura, Theoretical Background of the Fragment Molecular Orbital (FMO) Method and Its Implementation in GAMESS, in: D. Fedorov, K. Kitaura (Eds.), Fragment Molecular Orbital Method, CRC Press, Boca Raton, FL, 2009, pp. 5-36.
  * [41] D. Fedorov, K. Kitaura, Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method,  J. Phys. Chem. A 111 (2007) 6904-6914. https://doi.org/10.1021/jp0716740
  * [42] M.P. Mazanetz, E. Chudyk, D.G. Fedorov, Y. Alexeev, Applications of the Fragment Molecular Orbital Method to Drug Research, in: W. Zhang (Eds.), Computer-Aided Drug Discovery. Methods in Pharmacology and Toxicology. Humana Press, New York, NY, 2015.
  * [43] S. Eiler, M. Gangloff, S. Duclaud, D. Moras, M. Ruff, Overexpression, Purification, and Crystal Structure of Native ER$\alpha$LBD,  Protein Expression and Purification 22 (2001) 165-173. https://doi.org/10.1006/prep.2001.1409.
  * [44] D.A. Case, J.T. Berryman, R.M. Betz, D.S. Cerutti, T.E. Cheatham III, T.A. Darden, R.E. Duke, T.J. Giese, H. Gohlke, A.W. Goetz, N. Homeyer, S. Izadi, P. Janowski, J. Kaus, A. Kovalenko, T.S. Lee, S. LeGrand, P. Li, T. Luchko, R. Luo, B. Madej, K.M. Merz, G. Monard, P. Needham, H. Nguyen, H.T. Nguyen, I. Omelyan, A. Onufriev, D.R. Roe, A. Roitberg, R. Salomon-Ferrer, C.L. Simmerling, W. Smith, J. Swails, R.C. Walker, J. Wang, R.M. Wolf, X. Wu, D.M. York and P.A. Kollman (2015), AMBER 2015, University of California, San Francisco.
  * [45] J.M. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, Development and Testing of a General Amber Force Field,  J. Comput. Chem. 25 (2004) 1157-1174. https://doi.org/10.1002/jcc.20035.
  * [46] Gaussian 09, Revision A.1, M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, and D.J. Fox, Gaussian, Inc., Wallingford CT, 2009.
  * [47] C.I. Bayly, P. Cieplak, W. Cornell, P.A. Kollman, A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model,  J. Phys. Chem. 97 (1993) 10269-10280. https://doi.org/10.1021/j100142a004.
  * [48] H.G. Petersen, Accuracy and efficiency of the particle-mesh-ewald method,  J. Chem. Phys. 103 (1995) 3668-3679. https://doi.org/10.1063/1.470043.
  * [49] J-P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, Numerical integration of the cartesian equations of motion of a system with constraints molecular dynamics of n-alkanes,  J. Chem. Phys. 23 (1977) 327-341. https://doi.org/10.1016/0021-9991(77)90098-5.
  * [50] Ch.W. Hopkins, S. Le Grand, R.C. Walker, A.E. Roitberg, Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning,  J. Chem. Theory Comput. 11 (20015) 1864-1874. https://doi.org/10.1021/ct5010406.
  * [51] L.S.D. Caves, J.D. Evanseck, M. Karplus, Locally accessible conformations of proteins: Multiple molecular dynamics simulations of crambin,  Protein Science 7 (1998) 649-666. https://dx.doi.org/10.1002
  * [52] H.L. Gordon, R.L. Somorjai, Fuzzy Cluster Analysis of Molecular Dynamics Trajectories,  PROTEINS: Structure, Function, and Genetics 14 (1992) 249-264. https://doi.org/10.1002/prot.340140211.
  * [53] P.S. Shenkin, D.Q McDonald, Cluster analysis of molecular conformations,  Journal of Computational Chemistry 15 (1994) 899-916. https://doi.org/10.1002/jcc.540150811.
  * [54] M. Feig, J. Karanicolas, Ch. L. Brooks, III: MMTSB Tool Set, MMTSB NIH Research Resource,  The Scripps Research Institute (2001).
  * [55] W. Humphrey, A. Dalke, K. Schulten, VMD: Visual Molecular Dynamics,  J. Molec. Graphics 14 (1996) 33-38. https://doi.org/10.1016/0263-7855(96)00018-5.
  * [56] D.G. Fedorov, T. Ishida, K. Kitaura, Multilayer Formulation of the Fragment Molecular Orbital Method (FMO),  J. Phys. Chem. A 109 (2005) 2638-2646. https://doi.org/10.1021/jp047186z.
  * [57] I. Daidone, A. Amadei, Essential dynamics: foundation and applications,  WIREs Comput Mol Sci (2012) https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1099.
  * [58] A. Amadei, A.B.M. Linssen, H.J.C. Berendsen, Essential Dynamics of Proteins,  PROTEINS: Structure, Function, and Genetics 17 (1993) 412-425. https://doi.org/10.1002/prot.340170408.
  * [59] D.R. Roe, T.E. Cheatham, PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data,  J. Chem. Theory Comput. 9 (2013) 3084-3095. https://doi.org/10.1021/ct400341p.
  * [60] K. Fukuzawa, Y. Mochizuki, S. Tanaka, K. Kitaura, T. Nakano, Molecular Interactions between Estrogen Receptor and Its Ligand Studied by the ab Initio Fragment Molecular Orbital Method,  J. Phys. Chem. B 110 (2006) 16102-16110. https://doi.org/10.1021/jp065705n.
  * [61] A.H. Lima Costa, W.S. Clemente Jr, K.S. Bezerra, J.X. Lima Neto, E.L. Albuquerque, U.L. Fulco, Computational biochemical investigation of the binding energy interactions between an estrogen receptor and its agonists,  New J. Chem. 42 (2018) 19801-19810. https://doi.org/10.1039/c8nj03521k.