# Decoupling of dipolar and hydrophobic motions in biological membranes

Hanne S. Antila  Department of Theory and Bio-Systems, Max Planck Institute of
Colloids and Interfaces, 14424 Potsdam, Germany  Anika Wurl  Institute for
Physics, Martin-Luther University Halle-Wittenberg, 06120 Halle (Saale),
Germany  O. H. Samuli Ollila  Institute of Biotechnology, University of
Helsinki, 00014 Helsinki, Finland  Markus S. Miettinen  Department of Theory
and Bio-Systems, Max Planck Institute of Colloids and Interfaces, 14424
Potsdam, Germany  Tiago M. Ferreira<EMAIL_ADDRESS>Institute for Physics, Martin-Luther University Halle-Wittenberg

###### Abstract

Cells use homeostatic mechanisms to maintain an optimal composition of
distinct types of phospholipids in cellular membranes. The hydrophilic dipolar
layer at the membrane interface, composed of phospholipid headgroups,
regulates the interactions between cell membranes and incoming molecules,
nanoparticles, and viruses. On the other hand, the membrane hydrophobic core
determines membrane thickness and forms an environment for membrane-bound
molecules such as transmembrane proteins. A fundamental open question is to
what extent the motions of these regions are coupled and, consequently, how
strongly the interactions of lipid headgroups with other molecules depend on
the properties and composition of the membrane hydrophobic core. We combine
advanced solid-state nuclear magnetic resonance spectroscopy methodology with
high-fidelity molecular dynamics simulations to demonstrate how the rotational
dynamics of choline headgroups remain nearly unchanged (slightly faster) with
incorporation of cholesterol into a phospholipid membrane, contrasting the
well known extreme slowdown of the other phospholipid segments. Notably, our
results suggest a new paradigm where phospholipid headgroups interact as
quasi-freely rotating flexible dipoles at the interface, independent of the
properties in the hydrophobic region.

##  I Introduction

Out of the plethora of lipids found in nature, the most ubiquitous are
glycerophospholipids which consist of a glycerol backbone attached to two
hydrophobic fatty acid chains and a phosphodiester bridge connecting to a
hydrophilic headgroup Marsh (2013). Cells use vast amounts of energy, and rely
on complex synthetic pathways, to adjust and maintain the specific composition
of different types of phospholipid headgroups across cellular organelles
Harayama and Riezman (2018). This chemical homeostasis implies that the
headgroups play a key role in fundamental biological processes, and evidence
for phospholipid-specific functionality concerning compartmentalization,
signaling, transport, ion binding, peptide insertion, and regulation of
membrane protein function has been found Lemmon (2008); van Meer et al.
(2008); Harayama and Riezman (2018); Corradi et al. (2019). However, the
molecular details on how lipid composition of a cellular membrane connects to
its overall properties and to specific biological processes remain poorly
understood.

A key open question is to what degree the behavior of the water-facing
headgroups in biological membranes correlate with the properties of the acyl
chains in the membrane hydrophobic core Rajan et al. (1981); Huang and
Feigenson (1999); Roberts and Redfield (2004); Ali et al. (2007); Klauda et
al. (2008); Sivanandam et al. (2009); Alwarawrah et al. (2010); Leeb and
Maibaum (2018). Two limiting cases can be considered: 1) the conformational
ensemble and dynamics of the headgroup, though positionally connected through
the glycerol backbone, are uncoupled from the acyl chain region (freely
rotating/weak coupling limit), or 2) the orientation and dynamics of the
headgroup and hydrophobic acyl chains are strongly interdependent (strong
coupling limit). These two limiting scenarios will give rise to very distinct
biophysical behaviour. In the case of strong coupling, the headgroups in lipid
domains with ordered acyl chains and hindered dynamics, such as cholesterol-
induced lipid rafts Simons and Ikonen (1997, 2000), would exhibit slower
dynamics and possibly a different conformational ensemble. The acyl chain
behavior would then indirectly affect the interactions between lipid
headgroups and molecules in the aqueous media or within the membrane, such as
proteins or drugs. Such scenario is implicit, for example, in the popular
umbrella model for lipid–cholesterol interactions van Meer et al. (2008);
Huang and Feigenson (1999); Ali et al. (2007). In the weak coupling limit, the
behavior of lipid headgroups is similar irrespective of the acyl chain
structure, order, and dynamics, and consequently any cellular processes that
depend on the conformation and dynamics of the headgroups are unaffected by
the properties of the hydrophobic region.

We address the question of headgroup-tail decoupling using a
phosphatidylcholine (PC)–cholesterol bilayer system (the most abundant
phospholipid and sterol in eukaryotic cells)—a model cellular membrane from
which a wealth of both experimental and simulation data can be obtained.
Cholesterol is known to drive lateral heterogeneity and make membranes more
ordered (the so-called cholesterol condensing effect), which manifests as a
substantial increase in hydrophobic acyl chain C–H bond order parameters
($S_{\rm{CH}}$) in nuclear magnetic resonance (NMR) experiments Vist and Davis
(1990); Hofsäss et al. (2003); Tiburu et al. (2004); Warschawski and Devaux
(2005); Vermeer et al. (2007); Davis et al. (2009); Leftin and Brown (2011);
Ferreira et al. (2013a); Andersson et al. (2017). In contrast, the headgroup
and glycerol backbone order parameters are essentially unaffected up to the
highest cholesterol concentrations possible to incorporate in PC membranes
Brown and Seelig (1978); Ferreira et al. (2013b). The $\alpha$-carbon order
parameter of the choline headgroup, in particular, remains unchanged also upon
other bilayer perturbations that significantly affect acyl chain order
parameters, such as temperature, acyl chain composition, or membrane phase
(Table S1 in the supplementary information and references therein). On the
other hand, the choline headgroup orientation (and consequently the headgroup
order parameters) is highly sensitive to hydration level Ulrich and Watts
(1994), hydrostatic pressure Bonev and Morrow (1995), and the inclusion of
charges Altenbach and Seelig (1984); Scherer and Seelig (1989); Macdonald et
al. (1991a) or molecular dipoles in the membrane Bechinger and Seelig (1991a).
Therefore, a picture of a rotationaly decoupled headgroup, whose orientation
is independent of the hydrophobic region but can be affected by the
environment, emerges.

Although the C–H bond order parameters contain accurate information on
conformational ensembles, they do not convey how fast that ensemble is sampled
(conformational dynamics). The motional time-scales have been dominantly
assessed through spin-lattice relaxation ($R_{1}$) and spin-lattice relaxation
in the rotating frame ($R_{1\rho}$) NMR measurements Brown et al. (1979,
1983); Morrison and Bloom (1994); Klauda et al. (2008); Ferreira et al.
(2015); Sivanandam et al. (2009); Roberts M F (2009); Le Guernevé and Auger
(1995) which are sensitive to different (limited) time-scales depending on
experimental conditions and from which a physically meaningful change in the
dynamics (speedup vs. slowdown) can be challenging to interpret without
multiple measurements under different magnetic fields. Such measurements
demonstrate that the cholesterol-induced order in acyl chains is accompanied
by slower rotational dynamics not only of the acyl chains but also of the
glycerol backbone segments for which there is only a marginal conformational
change Sivanandam et al. (2009); Roberts M F (2009).

However, in a crucial contrast, the impact of cholesterol on the headgroup
motional time-scales, potentially occurring either via direct interaction or
through the observed glycerol backbone slowdown, has remained unclear. An
observation of increase of headgroup 13C $R_{1\rho}$ rates of DMPC upon
incorporation of cholesterol Le Guernevé and Auger (1995) suggests that the
cholesterol-induced slowdown of tail dynamics propagates to the phospholipid
headgroups although neither the statistical significance nor the quantitative
interpretation of the $R_{1\rho}$ increase in terms of physically meaningful
correlation times was provided. In stark contrast, comparison of 13C cross
polarization (CP) and refocused insensitive nuclear enhanced polarization
transfer (rINEPT) intensities suggest that the headgroup motional time-scales
remain unchanged even by the addition of 50% cholesterol Ferreira et al.
(2013a).

Here, we show that the dynamics of the PC headgroup are unaffected by
cholesterol, and consequently, that the motion of phospholipid headgroups is
decoupled from the hydrophobic region (the freely rotating limit). To this
end, we employ our novel NMR methodology Ferreira et al. (2015) where
segmental effective correlation times ($\tau_{\rm{e}}$) are determined from
solid-state NMR measurements of $R_{1}$, $R_{1\rho}$ and $S_{\rm{CH}}$. The
analysis of $\tau_{\rm{e}}$ values enables us to interpret the relaxation
rates in terms of a single, physically meaningful average time-scale for each
carbon where lower $\tau_{\rm{e}}$) value denotes slowdown and vise-versa.
Additionally, we decipher the origin of the decoupled motion by analysing
distinct all-atom (CHARMM36 and Slipids) and united-atom (Berger) MD
simulations which provide either realistic uncoupled (CHARMM36 and Slipids) or
non-realistic coupled (Berger) headgroup motions. The MD simulation models
indicate that the decoupled motion originates from dihedral rotations that are
present in all other glycerophospholipid types in addition to PCs. This
suggests that biological membranes have independent rotational dynamics of the
headgroups from the hydrophobic region, a feature that may be relevant in the
machinery of biological cell membrane processes.

##  II Materials and Methods

###  II.1 Sample Preparation

1-palmitoyl,2-oleoyl-$sn$-glycero-3-phosphocholine (POPC), cholesterol and
chloroform were purchased from Sigma-Aldrich. The samples were prepared by
mixing the lipids with chloroform and rapidly evaporating the organic solvent
under a nitrogen gas flow to obtain a homogeneous lipid film. Subsequently,
the lipid film was dried under vacuum overnight. The film was then hydrated in
a 0.5 ml tube by adding 50 %wt of water and manually mixing with a thin metal
rod multiple times alternated by sample centrifugation until a homogeneous
mixture was attained. The resulting mixture was then centrifuged into a KEL-F
Bruker insert with a sample volume of approximately 25 $\mu$l specifically
designed for solid-state NMR 4 mm rotors and left to equilibrate for at least
24 hours at room temperature before measurements.

###  II.2 NMR Experiments

The solid-state NMR experiments to measure $R_{1}$ and $R_{1\rho}$ were
performed on a Bruker Avance II-500 NMR spectrometer operating at a 13C Larmor
frequency of 125.78 MHz equipped with an E-free CP-MAS 4 mm (13C/31P/1H). The
R-PDLF measurements were performed on a Bruker Avance III 400 spectrometer
operating at a 1H Larmor frequency of 400.03 MHz equipped with a standard 4 mm
CP-MAS HXY probe. All experiments were performed under magic-angle spinning
(MAS) conditions at a rate of 5 kHz. The R-PDLF, $R_{1}$, and $R_{1\rho}$
experiments were performed as previously described in references Löser et al.
(2018); Ferreira et al. (2015). More details on experimental set up are given
in the supplementary information (SI).

###  II.3 MD Simulations

We performed MD simulations for two systems: a pure POPC bilayer and a bilayer
containing additional 50% cholesterol, both accompanied by enough water
molecules per lipid to result in fully hydrated bilayers. We used three lipid
MD models (force fields): CHARMM36 Klauda et al. (2010), Slipids Jämbeck and
Lyubartsev (2012), and Berger Ollila et al. (2007)/Höltje Höltje et al.
(2001); Ferreira et al. (2013c) together with either TIP3P Jorgensen et al.
(1983); MacKerell et al. (1998) or SPC Berendsen et al. (1981) water. The
choice of force fields was based on previous works where their ability to
capture structure Botan et al. (2015) and dynamics Antila et al. (2021) of
lipid headgroup and glycerol backbone was assessed against NMR measurables.
The simulations were performed using the GPU-version of Gromacs2020 Abraham et
al. (2015) MD engine, with sampling rate of 10 ps and maintaining 303 K
temperature. The list of all simulated systems, along with the trajectory
lengths and links to the freely available simulation data, is given in Table
1. Further details of the simulations are presented in the SI.

Table 1: Summary of simulated systems: the force fields used, numbers of POPC,
cholesterol, and water molecules, trajectory lengths, and access links to the
simulation files.

Force-field POPC/water+cholesterol | POPC/chol | water | length (ns) | files
---|---|---|---|---
CHARMM36 Klauda et al. (2010)/TIP3P MacKerell et al. (1998) | 122/0 | 4480 | 840 |  [51]
Slipids Jämbeck and Lyubartsev (2012)/TIP3P Jorgensen et al. (1983) | 122/0 | 4480 | 1200 |  [52]
Berger-POPC-07Ollila et al. (2007)/SPC Berendsen et al. (1981) | 256/0 | 10342 | 1200 |  [53]
CHARMM36 Klauda et al. (2010)/TIP3P MacKerell et al. (1998)+CHARMM36 Lim et al. (2012) | 122/122 | 9760 | 1240 |  [51]
Slipids Jämbeck and Lyubartsev (2012)/TIP3P Jorgensen et al. (1983)+Slipids Jämbeck and Lyubartsev (2013) | 122/122 | 9760 | 1200 |  [56]
Berger-POPC-07 Ollila et al. (2007)/SPC Berendsen et al. (1981)+Höltje-CHOL-13 Höltje et al. (2001); Ferreira et al. (2013c) | 256/256 | 20480 | 1200 |  [53]

Figure 1: Effect of cholesterol on the dynamics (panel B) and structure (panel
C) of headgroup ($\alpha$, $\beta$ and $\gamma$) and glycerol backbone
($g_{1}$, $g_{2}$, and $g_{3}$) carbons in POPC lipid membranes. (A) Chemical
structure of POPC with carbon labels and refocused-INEPT spectra from POPC
membranes with (red) and without (black) cholesterol. (B) 13C spin-lattice
relaxation rates, $R_{1}$, and spin-lattice relaxation in the rotating frame
rates, $R_{1\rho}$, showing the independent motion of the headgroup and
slowdown of the glycerol backbone upon cholesterol incorporation. A Larmor
frequency of 500 MHz for 1H nuclei and a spin lock field equal to 50 kHz were
used. The corresponding experimental decays for each data value are shown in
Figure S1. (C) Dipolar recoupling profiles acquired with R-PDLF spectroscopy
from POPC membranes with (red) and without (black) cholesterol. Similar
magnitudes of the splittings indicate structural independence of both the
headgroup and glycerol backbone on cholesterol incorporation. Note that the
line shapes are highly sensitive to the experimental setup and that the
relevant information on conformations is in the splittings, which are
proportional to $|S_{\rm{CH}}|$. (D) Overlayed lipid conformations in MD
simulations with (right) and without (left) cholesterol (yellow) illustrating
the experimental observations.

##  III Results

##  IV Experimental demonstration of the uncoupled motion

Figure 1 shows the effect of cholesterol on both the dynamics and structure of
POPC headgroup and glycerol backbone. Chemical shift resolution for all the
distinct carbons in the refocused-INEPT spectra displayed in Figure 1A is
enabled by simultaneous magic-angle spinning and heteronuclear decoupling. The
effect of cholesterol on the phospholipid dynamics is assessed by measuring
the $R_{1}$ and $R_{1\rho}$ values from POPC and POPC/cholesterol (1:1) multi-
lamellar vesicles and the effect on phospholipid structure by R-PDLF
spectroscopy. The relaxation rates for the headgroup and glycerol backbone are
shown in Figure 1B. The complete set of decays and relaxation rates measured
for all the phospholipid segments resolved in the 13C spectrum is given in
supplementary information Figures S1 (headgroup and glycerol backbone) and S2
(acyl chains). The dipolar splittings used to calculate the $S_{\rm{CH}}$
order parameters of headgroup and glycerol backbone are presented in Figure
1C. The resulting order parameters confirm the previously reported values
Ferreira et al. (2013a).

$R_{1}$ rates remain constant for both the glycerol backbone and the
headgroup, showing that the C–H bond motions with time-scales close to ns are
not affected by cholesterol for these carbons. On the other hand, the
$R_{1\rho}$ rates in the glycerol backbone (carbons $g_{1}$, $g_{2}$, and
$g_{3}$) increase by approximately a factor of two upon cholesterol addition.
In sharp contrast, the $R_{1\rho}$ values for the choline headgroup ($\alpha$,
$\beta$, and $\gamma$ segments) are unaffected by cholesterol and
significantly lower than in the glycerol backbone. The invariance of both the
headgroup carbon dipolar couplings (Figure 1C) and the relaxation rates
(Figure 1B) upon cholesterol incorporation implies that the conformational
ensemble and the time required to span all the available conformations is the
same irrespective of the glycerol backbone slowdown induced by cholesterol.
For acyl chains both structural ($S_{\rm{CH}}$) and dynamic observables
($R_{1}$ and $R_{1\rho}$) vary with incorporation of cholesterol (SI Figures
S2 and S3), signaling the expected ordering and slowdown in line previously
reported results Ferreira et al. (2013a); Sivanandam et al. (2009).

Figure 2: Impact of cholesterol (hollow bars: pure POPC, filled bars: 50%
POPC+50% cholesterol) on the effective correlation times, $\tau_{\rm{e}}$, of
different carbons in the headgroup and glycerol backbone of POPC quantified
experimentally and from lipid bilayer MD simulations with the CHARMM36,
Slipids and Berger force-fields. Note the different $y$-scales used on the
left and right plots to appreciate the significant difference of effective
correlation times for the choline headgroup (0.1-0.5 ns) and glycerol backbone
segments (2-5 ns).

To give an intuitive measure, sensitive to all time-scales, for the headgroup
and glycerol backbone dynamics, we quantified the effective correlation times
$\tau_{e}$. (Figure 2). To this end, we used the experimental data presented
in panels B and C of Figure 1 to calculate Ferreira et al. (2015)

$\tau_{e}=\frac{5R_{1\rho}-3.82R_{1}}{4\pi^{2}d_{\rm{CH}}^{2}N(1-S^{2}_{\rm{CH}})},$
(1)

where $N$ denotes the number of protons covalently bound to the carbon and the
coupling constant $d_{\rm{CH}}$ is approximately -22 kHz. The $\tau_{e}$
values of the choline headgroup are within 0.1–0.5 ns and remain constant
within the experimental accuracy upon addition of 50% cholesterol, while the
$\tau_{\rm{e}}$ values for the glycerol backbone slowdown with a factor of
approximately two, and are one order of magnitude slower than in the
headgroup.

##  V Comparison of experiments with MD simulations

Figure 3: Effect of cholesterol (hollow symbols: pure POPC, filled symbols:
50% POPC+50% cholesterol) on the POPC headgroup and glycerol backbone C–H bond
order parameters. The experimental values were determined by R-PDLF
spectroscopy. The grey areas show the range of $S_{\rm{CH}}$ values for PC
bilayer systems with and without cholesterol reported until date (see e.g.
References Leftin and Brown (2011); Ferreira et al. (2013a)). For comparing
with the effect on the acyl chains see Figure S3.  Figure 4: Effect of
cholesterol on internal structure and dynamics of POPC. (A) Dihedral angle
distributions for pure POPC membranes (darker thick lines) and
POPC/cholesterol membranes (lighter thin lines) from MD simulations using the
CHARMM36 (red), Slipids (blue) and Berger (green) force-fields. (B) Reduced
and normalized dihedral torsion autocorrelation functions (see Eq. 3) showing
their corresponding dihedral effective correlation times ($\tau_{c}^{\phi}$).
Note that the extremely short correlation times for the outermost right column
are due to the very narrow angle range accessible for this torsion dihedral.

Also included in Figure 2 are the $\tau_{\rm{e}}$ values calculated from three
sets of MD simulations using the CHARMM36, Slipids, and Berger force-fields.
Notably, CHARMM36 simulations reproduce the experimental $\tau_{\rm{e}}$
values of the choline headgroup, both with and without cholesterol, flawlessly
within experimental uncertainty. CHARMM36 also captures well the experimental
choline C-H bond order parameters which remain the same with the presence of
the sterol (Figure 3). For the glycerol backbone, CHARMM36 simulations
slightly overestimate the slowdown of the effective correlation times, but
give the best structural model among the three force-fields used. Slipids
simulations show the best agreement with experiments for the effect of
cholesterol on the glycerol backbone $\tau_{\rm{e}}$ values, although they
fail to capture the glycerol backbone structure, and overestimate the
effective correlation time for the $\gamma$ carbon.

The Berger force-field clearly produces the least realistic dynamics, giving a
significant overestimation of $\tau_{\rm{e}}$ for the choline headgroup
segments (both with and without cholesterol) and predicts an erroneous, large
(approximately 5-fold) cholesterol-induced slowdown of the choline headgroup
dynamics (see Figure S4 for $\tau_{e}$s from the Berger simulations drawn to
scale). Both the structural and dynamical force field properties observed here
are in line with the previously reported Botan et al. (2015); Antila et al.
(2021).

In contrast to the Berger model, both the CHARMM36 and Slipids models,
although not perfect, capture the key experimental observation in this work:
The structure and dynamics of the choline headgroup are not affected by
incorporation of cholesterol despite the increased acyl chain order and
hindered dynamics in the acyl chains and glycerol backbone, i.e, the headgroup
is uncoupled from the glycerol backbone in these models.

##  VI Effect of cholesterol on the time-scale of internal motions of the
headgroup and glycerol backbone

Knowing their ability to reproduce the NMR measurables, we proceed to exploit
the temporal and spatial resolution in the distinct MD models to gain insight
on the rotational dynamics of specific sites on the molecules as well as the
origin and degree of (un)coupling between the headgroup and the rest of the
phospholipid. To this end, Figure 4 shows the distributions of selected
headgroup and glycerol dihedral angles, $\phi$, and the corresponding dihedral
effective correlation times, $\tau_{c}^{\phi}$, which are extracted from
autocorrelation functions

$G^{\phi}(\tau)=\langle\cos\phi(t)\cos\phi(t+\tau)\rangle,$ (2)

where the angular brackets denote an average over time and over the number of
molecules in the system. We define the dihedral effective correlation time
($\tau_{c}^{\phi}$) as simply the area under the reduced normalised
autocorrelation function

$g^{\phi}(\tau)=\frac{G^{\phi}(\tau)/G^{\phi}(0)-\langle\cos\phi\rangle^{2}}{1-\langle\cos\phi\rangle^{2}},$
(3)

where $\langle\cos\phi\rangle^{2}$ is the value of $G^{\phi}(\tau)$ at
infinitely long $\tau$. $\tau_{c}^{\phi}$ provides a measure of how much time
is needed for dihedral motion to sample its angle distribution.

For the more realistic force-fields (CHARMM36 and Slipids) cholesterol does
not affect the dihedral distributions of POPC as expected from the
$S_{\rm{CH}}$ data alone. More interestingly, for both of these force-fields,
the dihedral angles over the phosphate linkage between the glycerol and the
alpha carbon, ($g_{3}$)O–P(O) and (O)P–O($\alpha$), and the glycerol backbone
dihedral $g_{2}$–$g_{3}$, have no angles which are unavailable, in contrast to
the other dihedrals analysed. Furthermore, the ($g_{3}$)O–P(O) and
(O)P–O($\alpha$) exhibit fast correlation times ($\tau_{c}^{\phi}\leq$0.5 ns)
which are very close to the effective correlation times in Figure 2.

The most notable differences between the more realistic simulations and Berger
simulations are the slower $\tau_{c}^{\phi}$ in Berger, and how these are
affected by cholesterol. While in CHARMM36 and Slipids only minor changes are
observed, with a slight speedup of the sampling by cholesterol, in the Berger
model these internal dynamics slowdown considerably.

##  VII Polar and azimuthal motion of the choline dipole and of the glycerol
backbone

To investigate the correlation between motions of headgroup and other parts of
lipid molecules, we quantified the autocorrelation functions of the polar and
azimuthal angles, $\theta$ and $\varphi$ (coordinate system where the
z-direction coincides with membrane normal), for a number of selected vectors
between intramolecular atomic pairs. The definition of the autocorrelation
functions $G^{\theta}(\tau)$ and $G^{\varphi}(\tau)$ are the same as in Eq. 2
but using $\theta$ and $\varphi$ as angles, respectively. Note that for
$G^{\theta}(\tau)$ a non-zero plateau at the long $\tau$ is expected since the
different $\theta$ angles are not equally likely to occur. On the other hand,
$G^{\varphi}(\tau)$ is always zero at long $\tau$ due to the lipid uniaxial
motion.

In Figure 5 we show these correlation functions for the choline dipole
orientation, P$\rightarrow$N, and for the interatomic vector connecting the
carbonyl carbons in the sn-1 and sn-2 positions, calculated from the most
realistic force-field (CHARMM36). The $\varphi$ autocorrelation functions
clearly show the contrasting effects of cholesterol on these vectors. While
for the choline dipole a speedup of reorientational motion is observed, the
$\varphi$ dynamics of the vector connecting the carbonyl carbons become slower
by almost an order of magnitude. To extract effective correlation times
$\tau_{c}^{\theta}$ and $\tau_{c}^{\varphi}$, we again integrate the reduced
and normalised autocorrelation functions. While cholesterol induces more than
a 4-fold slowdown for the $\varphi$ dynamics of the vector connecting carbonyl
carbons, the correlation times of the P$\rightarrow$N orientation remain
essentially the same. The complete set of reduced autocorrelations functions
analysed is given in SI Figures S7-S13 together with $\theta$ distributions,
and $\tau_{c}^{\theta}$ and $\tau_{c}^{\varphi}$ values.

Figure 5: The effect of cholesterol on the time-scales for the reorientations
of P$\rightarrow$N and (sn-1)O=C$\rightarrow$C=O(sn-2) vectors over the
spherical angles $\varphi$ and $\theta$. The autocorrelation functions shown
here are described in Equation 2.

##  VIII Discussion

Our experimental and MD simulation results show that the phospholipid
headgroup conformational ensemble and dynamics remain unaffected by addition
of 50% cholesterol to the lipid membrane, despite the significant acyl chain
ordering and reduction in both the acyl chain and glycerol backbone dynamics.
Therefore, our results do not support models that contain interdependence
between the structure or dynamics of the hydrophobic and hydrophilic regions
of cellular or model phospholipid membranes.

The observed slowdown of the glycerol backbone upon addition of cholesterol
(Figures 2 and S5) arises from the longer timescales to which $R_{1\rho}$ is
sensitive to. The internal motions of the glycerol backbone are not affected
by cholesterol in the most realistic MD simulations used (CHARMM36 and
Slipids, see Figure 4) in line with the invariance of the glycerol backbone
$R_{1}$ values (Figure 1). Therefore, the slower glycerol backbone dynamics
induced by cholesterol most likely arises from a slowdown of the rotational
diffusion of the whole phospholipid body as previously suggested by Roberts et
al. Sivanandam et al. (2009); Roberts M F (2009), rather than restrictions in
internal dynamics.

The independence of headgroup and hydrophobic chain motions must result from a
set of fast internal rotations around phospholipid bonds with specific
orientations that decouple these motions. The MD simulations in best agreement
with the NMR experiments show a high flexibility for dihedral angles in the
headgroup region with a wide range of accessible conformations (Figures 4, S6
as well as Ref. Bacle et al. (2021)). From the set of dihedral distribution
functions, one clearly observes the highly flexible nature of the
($g_{3}$)O–P, P–O($\alpha$) and $g_{2}$-$g_{3}$ dihedral angles with all
angles over a complete dihedral rotation having non-zero probability in
contrast to the remaining torsions. This applies for all three force-fields,
though Berger has a less even distribution than CHARMM36 and Slipids. For
CHARMM36 and Slipids, the correlation times for the rotations around the
($g_{3}$)O–P and P–O($\alpha$) bonds are lower than 0.5 ns (Figure 4) and very
close to the effective correlation times measured for the C–H bonds from the
$\alpha$ and $\beta$ carbons. These $\tau_{c}^{\phi}$ values slightly decrease
with the addition of cholesterol, i.e. cholesterol induces a slight speed-up
of the torsion dynamics for these particular dihedrals, most likely because
fewer steric hindrances are present due to an increased average distance
between headgroups.

The highly flexible dihedral rotations around the phosphate P–O bonds, as well
as the g2–g3 torsion, are much faster than the transverse and longitudinal
rotational diffusion of the molecular frame whose correlation times have been
approximated to be 10–20 ns and 100 ns, respectively Klauda et al. (2008). The
most probable $\theta$ angle of the the (g3–)O–P and g2–g3 bonds is less than
10∘ (Figures S8 and S11) which is very close to be parallel with the bilayer
normal axis. A flexible, fast-rotating dihedral aligned with the membrane
normal leads to decoupling of headgroup from the longitudinal rotational
diffusion of the molecular frame. Decoupling from the transverse rotational
diffusion is mostly enabled by the fast motion over the P–O(–Cα).

The torsions around the phosphate group and g2–g3 dihedrals act analogously to
a frictionless spherical-joint which decouples the choline headgroup structure
and dynamics from the glycerol backbone. This is in line with the previous
analysis Klauda et al. (2008) where a partial decoupling of the headgroup from
the main phospholipid body due to the rotation of a phosphate dihedral was
suggested based on a comparison of CHARMM C27r MD simulation of pure
dipalmitoylphosphatidylcholine bilayers to 31P-NMR R1 data under several
magnetic fields. Here, we demonstrate that such decoupling is strong enough to
prevent the propagation of the slowdown effect of cholesterol from the acyl
chains and glycerol backbone to the headgroup.

We base our molecular interpretation of the decoupled motion on the CHARMM36
force field, which gives the best overall description for the headgroup and
glycerol backbone structure and dynamics among the available models (Figures 2
and 3, and Refs. Botan et al. (2015); Antila et al. (2021)). However, not all
the NMR observables calculated are within experimental errors even in CHARMM36
simulations and we cannot exclude the possibility that improved future models
correctly capturing the decoupling effect, effective correlation times, and
structural order parameters may give an alternative molecular interpretation.

The headgroup decoupling is not observed in Berger force-field, although it
also provides fast dynamics over the phosphate group dihedrals (Figure 4).
This is most likely due to an overestimation of the attractive interaction
between cholesterol and the choline group. Such interaction has been
interpreted previously as a consequence of the so-called umbrella effect where
a reorientation of the headgroup due to presence of cholesterol is often
assumed Alwarawrah et al. (2012). However, the combination of MD simulations
and experiments presented here indicates that such interaction is artefactual
and that both the orientation and dynamics of the headgroup are unaffected by
the sterol presence. The implicit assumption in the umbrella model of a
cholesterol effect on headgroup reorientation, either through a change of the
conformational ensemble or a change of dynamics, is not supported by our
experimental results or the more realistic MD models.

Correlation functions of P$\rightarrow$N and other intramolecular vectors,
calculated from CHARMM36 simulations, further support the idea of decoupled
motions between headgroup and the acyl chains. Cholesterol induces a
significant slowdown of the reorientation of the interatomic vectors between
atoms belonging to the glycerol backbone and acyl chain segments (Figures
S7-S13), while the effect on the P$\rightarrow$N vector, representing the
choline dipole, is negligible with only a slight speed up of the dynamics most
likely due to the increase of the distance between phospholipids headgroups.

The molecular description suggested here has rather strong implications for
membrane biophysics and should motivate a number of additional experiments and
simulations. It implies that the dipolar surface of glycerophospholipid
bilayers consists of freely rotating dipoles with timescales faster than 2 ns
that do not depend on the dynamics of the acyl chains or glycerol backbone.
The time-scale of reorientation of the dipoles is expected to influence the
interaction of the headgroups with charged molecules, e.g. proteins, that
approach the lipid biomembrane. For instance, it is known that the tilt of the
headgroup dipole is highly sensitive to membrane surface charge Scherer and
Seelig (1989). Under a positive surface charge, the headgroups tilt to a more
upright orientation (increase of the alpha and beta $S_{\mathrm{CH}}$ values)
and vice-versa for a negative surface charge due to the charge-dipole
electrostatic interactions. The results presented here suggest that the
phosphate and the g2–g3 dihedrals enable an unconstrained response of
headgroups to the electrostatic field and effectively uncouple the
interactions occurring in the membrane surface from the hydrophobic region.
Although we only investigate here PC headgroups, it is foreseable that the
decoupling applies to all other glycerophospholipids since the same  molecular
bearings are present irrespectively of the substituent headgroup Marsh (2013).

In summary, our results suggest that for describing the dipolar interactions
at the surface of membranes, the hydrophobic structure may be neglected to a
good approximation and that the relevant headgroup physics lie on the
electrostatic interactions—which is remarkably useful considering the complex
molecular arrangement in the hydrophobic region of biological membranes.

##  IX Acknowledgments

O.H.S.O acknowledges CSC – IT Center for Science for computational resources
and Academy of Finland (grants 315596 and 319902) for financial support.
H.S.A. gratefully acknowledges financial support from the Osk. Huttunen
Foundation, Finnish Academy of Science and Letters (Foundations’ Post Doc
Pool), Instrumentarium Science Foundation, and the Alexander von Humboldt
Foundation. T.M.F. was supported by the Ministry of Economics, Science and
Digitalisation of the State of Saxony-Anhalt, Germany.

T.M.F. greatly acknowledges Kay Saalwächter and Alexey Krushelnitsky for
invaluable support and discussions.

##  X Competing interests

The authors declare no competing interests.

## References

  * Marsh (2013) D. Marsh, _Handbook of Lipid Bilayers_ (CRC press, New York, 2013).
  * Harayama and Riezman (2018) T. Harayama and H. Riezman, Nat. Rev. Mol. Cell Biol. 19, 281 (2018), URL https://doi.org/10.1038/nrm.2017.138.
  * Lemmon (2008) M. A. Lemmon, Nat. Rev. Mol. Cell Biol. 9, 99 (2008), URL https://doi.org/10.1038/nrm2328.
  * van Meer et al. (2008) G. van Meer, D. R. Voelker, and G. W. Feigenson, Nat. Rev. Mol. Cell Biol. 9, 112 (2008), URL http://dx.doi.org/10.1038/nrm2330.
  * Corradi et al. (2019) V. Corradi, B. I. Sejdiu, H. Mesa-Galloso, H. Abdizadeh, S. Y. Noskov, S. J. Marrink, and D. P. Tieleman, Chem. Rev. 119, 5775 (2019), URL https://doi.org/10.1021/acs.chemrev.8b00451.
  * Rajan et al. (1981) S. Rajan, S. Y. Kang, H. S. Gutowsky, and O. E, J Biol Chem. 256, 1160 (1981), URL http://www.jbc.org/content/256/3/1160.short.
  * Huang and Feigenson (1999) J. Huang and G. W. Feigenson, Biophys. J. 76, 2142 (1999), URL https://pubmed.ncbi.nlm.nih.gov/10096908.
  * Roberts and Redfield (2004) M. Roberts and A. Redfield, Proc. Natl. Aca. Sci. USA 101, 17066 (2004), URL https://www.pnas.org/content/101/49/17066.
  * Ali et al. (2007) M. Ali, K. Cheng, and J. Huang, Proc. Natl. Acad. Sci. USA 104, 5372 (2007), URL https://www.pnas.org/content/104/13/5372.
  * Klauda et al. (2008) J. B. Klauda, M. F. Roberts, A. G. Redfield, B. R. Brooks, and R. W. Pastor, Biophys J. 94, 3074 (2008), URL http://www.sciencedirect.com/science/article/pii/S0006349508704648.
  * Sivanandam et al. (2009) V. N. Sivanandam, J. Cai, A. G. Redfield, and M. F. Roberts, J. Am. Chem. Soc. 131, 3420 (2009), URL https://doi.org/10.1021/ja808431h.
  * Alwarawrah et al. (2010) M. Alwarawrah, J. Dai, and J. Huang, J. Phys. Chem. B 114, 7516 (2010), URL https://pubmed.ncbi.nlm.nih.gov/20469902.
  * Leeb and Maibaum (2018) F. Leeb and L. Maibaum, Biophys. J. 115, 2179 (2018), URL https://pubmed.ncbi.nlm.nih.gov/30447996.
  * Simons and Ikonen (1997) K. Simons and E. Ikonen, Nature 387, 569 (1997), URL https://doi.org/10.1038/42408.
  * Simons and Ikonen (2000) K. Simons and E. Ikonen, Science 290, 1721 (2000), URL https://science.sciencemag.org/content/290/5497/1721.
  * Vist and Davis (1990) M. R. Vist and J. H. Davis, Biochemistry 29, 451 (1990), URL https://doi.org/10.1021/bi00454a021.
  * Hofsäss et al. (2003) C. Hofsäss, E. Lindahl, and O. Edholm, Biophys. J. 84, 2192 (2003), URL http://www.sciencedirect.com/science/article/pii/S0006349503750255.
  * Tiburu et al. (2004) E. K. Tiburu, P. C. Dave, and G. A. Lorigan, Magn. Reson. Chem. 42, 132 (2004), URL https://onlinelibrary.wiley.com/doi/abs/10.1002/mrc.1324.
  * Warschawski and Devaux (2005) D. E. Warschawski and P. F. Devaux, Eur. Biophys. J. 34, 987 (2005), URL https://doi.org/10.1007/s00249-005-0482-z.
  * Vermeer et al. (2007) L. S. Vermeer, B. L. de Groot, V. Reat, A. Milon, and J. Czaplicki, Eur. Biophys. J. Biophy. 36, 919 (2007), URL https://doi.org/10.1007/s00249-007-0192-9.
  * Davis et al. (2009) J. H. Davis, J. J. Clair, and J. Juhasz, Biophys. J. 96, 521 (2009), URL https://doi.org/10.1016/j.bpj.2008.09.042.
  * Leftin and Brown (2011) A. Leftin and M. F. Brown, Biochim. Biophys. Acta-Biomembr. 1808, 818 (2011), URL https://doi.org/10.1016/j.bbamem.2010.11.027.
  * Ferreira et al. (2013a) T. M. Ferreira, F. Coreta-Gomes, O. H. S. Ollila, M. J. Moreno, W. L. C. Vaz, and D. Topgaard, Phys. Chem. Chem. Phys. 15, 1976 (2013a), URL http://dx.doi.org/10.1039/C2CP42738A.
  * Andersson et al. (2017) J. M. Andersson, C. Grey, M. Larsson, T. M. Ferreira, and E. Sparr, P. Natl. Acad. Sci. USA 114, E3592 (2017), URL https://doi.org/10.1073/pnas.1701239114.
  * Brown and Seelig (1978) M. Brown and J. Seelig, Biochemistry 17, 381 (1978), URL https://doi.org/10.1021/bi00595a029.
  * Ferreira et al. (2013b) T. M. Ferreira, D. Bernin, and D. Topgaard (Academic Press, 2013b), vol. 79 of _Annual Reports on NMR Spectroscopy_ , pp. 73 – 127, URL http://www.sciencedirect.com/science/article/pii/B9780124080980000033.
  * Ulrich and Watts (1994) A. S. Ulrich and A. Watts, Biophys. J. 66, 1441 (1994), URL https://doi.org/10.1016/S0006-3495(94)80934-8.
  * Bonev and Morrow (1995) B. B. Bonev and M. R. Morrow, Biophys J. 69, 518 (1995), URL https://doi.org/10.1016/S0006-3495(95)79925-8.
  * Altenbach and Seelig (1984) C. Altenbach and J. Seelig, Biochemistry 23, 3913 (1984), URL https://doi.org/10.1021/bi00312a019.
  * Scherer and Seelig (1989) P. G. Scherer and J. Seelig, Biochemistry 28, 7720 (1989).
  * Macdonald et al. (1991a) P. M. Macdonald, J. Leisen, and F. M. Marassi, Biochemistry 30, 3558 (1991a), URL http://dx.doi.org/10.1021/bi00228a029.
  * Bechinger and Seelig (1991a) B. Bechinger and J. Seelig, Biochemistry 30, 3923 (1991a), URL https://doi.org/10.1021/bi00230a017.
  * Brown et al. (1979) M. F. Brown, J. Seelig, and U. Häberlen, J. Chem. Phys. 70, 5045 (1979), URL https://doi.org/10.1063/1.437346.
  * Brown et al. (1983) M. F. Brown, A. A. Ribeiro, and G. D. Williams, Proc. Natl. Acad. Sci. USA 80, 4325 (1983), URL https://doi.org/10.1073/pnas.80.14.4325.
  * Morrison and Bloom (1994) C. Morrison and M. Bloom, J. Chem. Phys. 101, 749 (1994), URL https://doi.org/10.1063/1.468491.
  * Ferreira et al. (2015) T. M. Ferreira, O. H. S. Ollila, R. Pigliapochi, A. P. Dabkowska, and D. Topgaard, J. Chem. Phys. 142, 044905 (2015), URL https://doi.org/10.1063/1.4906274.
  * Roberts M F (2009) M. U. Roberts M F, Redfield A G, Biophys J. 97, 132 (2009), URL https://doi.org/10.1016/j.bpj.2009.03.057.
  * Le Guernevé and Auger (1995) C. Le Guernevé and M. Auger, Biophys J. 68, 1952 (1995), URL https://doi.org/10.1016/S0006-3495(95)80372-3.
  * Löser et al. (2018) L. Löser, K. Saalwächter, and T. Ferreira, Phys. Chem. Chem. Phys. 20, 9751 (2018), URL http://dx.doi.org/10.1039/C8CP01012A.
  * Klauda et al. (2010) J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O’Connor, D. J. Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. MacKerell Jr, and R. W. Pastor, J. Phys. Chem. B 114, 7830 (2010).
  * Jämbeck and Lyubartsev (2012) J. P. M. Jämbeck and A. P. Lyubartsev, J. Chem. Theory Comput. 8, 2938 (2012).
  * Ollila et al. (2007) S. Ollila, M. T. Hyvönen, and I. Vattulainen, J. Phys. Chem. B 111, 3139 (2007).
  * Höltje et al. (2001) M. Höltje, T. Förster, B. Brandt, T. Engels, W. von Rybinski, and H.-D. Höltje, Biochim. Biophys. Acta 1511, 156 (2001).
  * Ferreira et al. (2013c) T. M. Ferreira, F. Coreta-Gomes, O. H. S. Ollila, M. J. Moreno, W. L. C. Vaz, and D. Topgaard, Phys. Chem. Chem. Phys. 15, 1976 (2013c), URL http://dx.doi.org/10.1039/C2CP42738A.
  * Jorgensen et al. (1983) W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983).
  * MacKerell et al. (1998) A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, et al., J. Phys. Chem. B 102, 3586 (1998).
  * Berendsen et al. (1981) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, _Interaction Models for Water in Relation to Protein Hydration_ (Springer Netherlands, Dordrecht, 1981), pp. 331–342.
  * Botan et al. (2015) A. Botan, F. Favela-Rosales, P. F. J. Fuchs, M. Javanainen, M. Kanduč, W. Kulig, A. Lamberg, C. Loison, A. Lyubartsev, M. S. Miettinen, et al., J. Phys. Chem. B 119, 15075 (2015), URL https://doi.org/10.1021/acs.jpcb.5b04878.
  * Antila et al. (2021) H. S. Antila, T. M. Ferreira, O. H. S. Ollila, and M. S. Miettinen, Journal of Chemical Information and Modeling 61, 938 (2021), eprint https://doi.org/10.1021/acs.jcim.0c01299, URL https://doi.org/10.1021/acs.jcim.0c01299.
  * Abraham et al. (2015) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1, 19 (2015).
  * Antila (2020a) H. Antila, _CHARMM36 simulations of pure POPC bilayer and a bilayer containing 50% of cholesterol, T=303K_ (2020a), URL https://doi.org/10.5281/zenodo.4071264.
  * Antila (2020b) H. Antila, _Slipids simulation of pure POPC bilayer at T=303K_ (2020b), URL https://doi.org/10.5281/zenodo.4352287.
  * Antila (2021) H. Antila, _MD simulation of 50% POPC+50% cholesterol bilayer using Berger force field at T=303K_ (2021), URL https://doi.org/10.5281/zenodo.4537794.
  * Lim et al. (2012) J. B. Lim, B. Rogaski, and J. B. Klauda, J. Phys. Chem. B 116, 203 (2012).
  * Jämbeck and Lyubartsev (2013) J. P. M. Jämbeck and A. P. Lyubartsev, J. Chem. Theory Comput. 9, 774 (2013), eprint https://doi.org/10.1021/ct300777p, URL https://doi.org/10.1021/ct300777p.
  * Antila (2020c) H. Antila, _Slipids simulations of a bilayer containing 50% of cholesterol, T=303K_ (2020c), URL https://doi.org/10.5281/zenodo.4349753.
  * Bacle et al. (2021) A. Bacle, P. Buslaev, R. G. Fandiño, F. Favela-Rosales, T. M. Ferreira, P. F. Fuchs, I. Gushchin, M. Javanainen, A. M. Kiirikki, J. J. Madsen, et al., Submitted (2021), URL https://github.com/NMRLipids/NMRlipidsIVPEandPGmanuscript/blob/master/manuscriptPGPE.pdf.
  * Alwarawrah et al. (2012) M. Alwarawrah, J. Dai, and J. Huang, Journal of Chemical Theory and Computation 8, 749 (2012).
  * Ollila (2014) O. H. S. Ollila, _MD simulation trajectory and related files for POPC/cholesterol (50 mol%) bilayer (Berger model delivered by Tieleman, modified Höltje, Gromacs 4.5)_ (2014), URL http://dx.doi.org/10.5281/zenodo.13285.
  * Bacle and Fuchs (2018) A. Bacle and P. F. Fuchs, _Berger pure popc md simulation (300 k - 300ns - 1 bar)_ (2018), URL https://doi.org/10.5281/zenodo.1402417.
  * Berendsen et al. (1984) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, The Journal of Chemical Physics 81, 3684 (1984), eprint https://doi.org/10.1063/1.448118, URL https://doi.org/10.1063/1.448118.
  * Parrinello and Rahman (1981) M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
  * Nose (1984) S. Nose, Mol. Phys. 52, 255 (1984).
  * Hoover (1985) W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
  * Essman et al. (1995) U. L. Essman, M. L. Perera, M. L. Berkowitz, T. Larden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995).
  * Hess (2008) B. Hess, Journal of Chemical Theory and Computation 4, 116 (2008).
  * Miyamoto and Kollman (1992) S. Miyamoto and P. A. Kollman, J. Comput. Chem 13, 952 (1992).
  * Macdonald et al. (1991b) P. M. Macdonald, J. Leisen, and F. M. Marassi, Biochemistry 30, 3558 (1991b), URL https://doi.org/10.1021/bi00228a029.
  * Bechinger and Seelig (1991b) B. Bechinger and J. Seelig, Chemistry and Physics of Lipids 58, 1 (1991b), URL http://www.sciencedirect.com/science/article/pii/000930849190105K.
  * Gally et al. (1975) H. U. Gally, W. Niederberger, and J. Seelig, Biochemistry 14, 3647 (1975), URL https://doi.org/10.1021/bi00687a021.

##  XI Supplementary Information

##  XII Solid-state NMR experiments

###  XII.1 R-PDLF experiments

A total of 32 points in the indirect dimension with increments equal to two
R18 blocks; SPINAL64 was used for proton decoupling during 13C acquisition,
with a nutation frequency of approximately 50 kHz, a total acquisition time of
0.07 s and a spectral width of 200 ppm; the rINEPT pulses were set at a
nutation frequency of 78.12 kHz.

###  XII.2 $R_{1}$ and $R_{1\rho}$ experiments

RF $\pi/2$ and $\pi$ pulses were set to a nutation frequency of 63.45 kHz.
TPPM was used for proton decoupling during 13C acquisition, with a nutation
frequency of approximately 50 kHz, a total acquisition time of 0.1 s, recycle
delay of 10 s and a spectral width of 140 ppm. The spin-lock frequency for
$R_{1\rho}$ was 50 kHz.

For quantifying $R_{1}$ and $R_{1\rho}$ for a given carbon segment, we
determined the decay over the indirect dimension by fitting gaussian
lineshapes in the direct dimension and using the analytic areas of the fitted
functions. The decay was then fitted with a single exponential decay and the
error bounds for both the $R_{1}$ and $R_{1\rho}$ values presented are the 95
% confidence bounds from these fits.

##  XIII MD simulations

The force field parameters were acquired from CHARMM-GUI (CHARMM36), the
Slipids web page (http://www.fos.su.se/~sasha/SLipids/, the pre 2020 version),
and from existing simulations in zenodo repository Ollila (2014); Bacle and
Fuchs (2018) (Berger). Similarly, the initial configurations, where the
bilayers spanned the x–y plane of a rectangular simulation box with periodic
boundary condition, were obtained either from CHARMM-GUI or constructed from
existing trajectories Ollila (2014); Bacle and Fuchs (2018). The initial
configurations were then relaxed using the steepest decent algorithm (5000
steps), followed by 2.5 ns + 41 ns of NPT simulation with semi-isotropic
pressure coupling maintaining 1 bar pressure. In the first 2.5 ns the pressure
set using the Berendsen barostat Berendsen et al. (1984) with coupling
constant 1 ps while the rest of the 41 ns pre-equilibration was done by
utilizing the Parrinello–Rahman barostat Parrinello and Rahman (1981) with 1
ps coupling constant. The latter barostat settings where then used for the
following 0.84 $\mu$s-1.24 $\mu$s of NPT production runs used for analysis
(See Table 1 in main text). The temperature in the simulations was set to 303
K using the Nose-Hoover thermostat Nose (1984); Hoover (1985) with coupling
constant 1 ps. The electrostatic interactions were modelled utilizing the
particle-mesh-Ewald algorithm Essman et al. (1995). For CHARMM36 and Slipids
simulations a cutoff of 1.2 nm was used to separate the real and reciprocal
space contributions to electrostatics whereas the Lennard-Jones potentials
were smoothly fixed to zero between cutoffs 1 nm and 1.2 nm. In Berger
simulations we utilized verlet list scheme with cutoff 1 nm for both
electrostatic and Lennard-Jones interactions. The leap-frog algorithm with
time step 0.002 ps was used to integrate the movement of the particles. To
avoid the need for a shorter time-step, the fast vibrations of the C–H bonds
of the lipids where removed using the LINCS algorithm Hess (2008) and SETTLE
Miyamoto and Kollman (1992) constrains were used for water. The resulting
simulation trajectories were analysed using in-house scripts. The exact
methodology utilized for extracting the NMR measurables, and their error
estimates, from the simulation is detailed in Ref. 49.

System | method | conditions |  $T$ / ∘C | Phase | $S^{\alpha}_{\rm{CH}}$ | $-S^{\beta}_{\rm{CH}}$ | ref.
---|---|---|---|---|---|---|---
DMPC |  2H NMR | 150 mM NaCl, 10 mM HEPES | 30 | Lα | 0.048 | 0.046 | Macdonald et al. (1991b)
DMPC |  2H NMR | 150 mM NaCl, 10 mM HEPES | 15 | Lβ | 0.051 | 0.062 | Macdonald et al. (1991b)
DOPC |  2H NMR | water only | 30 | Lα | 0.049 | 0.029 | Ulrich and Watts (1994)
POPC |  2H NMR | NaCl saturated solution | 23 | Lα | 0.048 | 0.044 | Bechinger and Seelig (1991b)
DPPC |  2H NMR | water only | 45-90 | Lα | 0.048 | 0.046-0.024 | Gally et al. (1975)
DPPC/chol. (1:1) |  2H NMR | 0.2 M sodium acetate/acetic acid | 10-70 | Lα | 0.048/0.048-0.040 | 0.03-0.024 | Brown and Seelig (1978)
POPC/chol. (1:1) | R-PDLF | water only | 30 | Lα | 0.052 | 0.034 | here
POPC | R-PDLF | water only | 30 | Lα | 0.052 | 0.040 | here
Table S1: Previously published $\alpha$ and $\beta$ C–H bond order parameters
from 2H NMR spectroscopy for a number of phosphatidylcholine lamellar systems
together with values reported here using 1H-13C dipolar recoupling on POPC and
POPC/cholesterol (1:1).  Figure S1: The $R_{1}$ and $R_{1\rho}$ decays
measured for the headgroup and glycerol backbone carbons in the POPC (blue)
and POPC/cholesterol (red) systems. Each point corresponds to the integral
determined from a Gaussian fit of the corresponding 13C peak in the high
resolution chemical shift spectrum acquired under MAS of 5 kHz. The spin lock
field for the $R_{1\rho}$ measurement was 50 kHz and the 13C Larmor frequency
was 125 MHz. Figure S2: The $R_{1}$ and $R_{1\rho}$ decays for the acyl chain
carbons measured in the POPC (blue) and POPC/cholesterol (red) systems. Each
point corresponds to the integral determined from a gaussian fit of the
corresponding 13C peak in the high resolution chemical shift spectrum acquired
under MAS of 5 kHz. The spin lock field for the $R_{1\rho}$ measurement was 50
kHz and the 13C Larmor frequency was 125 MHz. Figure S3: Effect of cholesterol
on the C–H bond order parameter magnitudes, $|S_{\rm{CH}}|$, of different
segments in the acyl chains of POPC measured experimentally and calculated
from lipid bilayer MD simulations with the CHARMM36 force-field. Figure S4:
Same as Figure 2 but with proper scale for visualizing the Berger data.
Figure S5: Effect of cholesterol on the effective correlation times,
$\tau_{\rm{e}}$, of different segments in the acyl chains of POPC calculated
experimentally and from lipid bilayer MD simulations with the CHARMM36 force
field.  Figure S6: Dihedral angle distributions from the POPC (thick black
lines) and POPC/cholesterol (thin colored lines) MD simulations using the
CHARMM36 (top), Slipids (middle) and Berger (bottom) force-fields.  Figure S7:
Orientation distribution, $p(\theta)$ (outermost left), and auto-correlation
functions for the polar, $\theta$ (middle row), and azimuthal, $\varphi$
(outermost right), angles of vector P$\rightarrow$N from POPC (top) and
POPC/cholesterol (bottom) CHARMM36 MD simulations. The grey lines represent
data for each individual phospholipid molecule in the simulation. $\theta$ and
$\varphi$ are the spherical angles in the laboratory coordinate frame defined
by the simulation box axes.  Figure S8: Orientation distribution, $p(\theta)$
(outermost left), and auto-correlation functions for the polar, $\theta$
(middle row), and azimuthal, $\varphi$ (outermost right), angles of vector
($g_{3}$-)O$\rightarrow$P from POPC (top) and POPC/cholesterol (bottom)
CHARMM36 MD simulations. The grey lines represent data for each individual
phospholipid molecule in the simulation. $\theta$ and $\varphi$ are the
spherical angles in the laboratory coordinate frame defined by the simulation
box axes.  Figure S9: Orientation distribution, $p(\theta)$ (outermost left),
and auto-correlation functions for the polar, $\theta$ (middle row), and
azimuthal, $\varphi$ (outermost right), angles of vector
($g_{1}$-)O$\rightarrow$O(-$g_{2})$ from POPC (top) and POPC/cholesterol
(bottom) CHARMM36 MD simulations. The grey lines represent data for each
individual phospholipid molecule in the simulation. $\theta$ and $\varphi$ are
the spherical angles in the laboratory coordinate frame defined by the
simulation box axes.  Figure S10: Orientation distribution, $p(\theta)$
(outermost left), and auto-correlation functions for the polar, $\theta$
(middle row), and azimuthal, $\varphi$ (outermost right), angles of vector
$g_{1}\rightarrow g_{3}$ from POPC (top) and POPC/cholesterol (bottom)
CHARMM36 MD simulations. The grey lines represent data for each individual
phospholipid molecule in the simulation. $\theta$ and $\varphi$ are the
spherical angles in the laboratory coordinate frame defined by the simulation
box axes.  Figure S11: Orientation distribution, $p(\theta)$ (outermost left),
and auto-correlation functions for the polar, $\theta$ (middle row), and
azimuthal, $\varphi$ (outermost right), angles of vector $g_{2}\rightarrow
g_{3}$ from POPC (top) and POPC/cholesterol (bottom) CHARMM36 MD simulations.
The grey lines represent data for each individual phospholipid molecule in the
simulation. $\theta$ and $\varphi$ are the spherical angles in the laboratory
coordinate frame defined by the simulation box axes.  Figure S12: Orientation
distribution, $p(\theta)$ (outermost left), and auto-correlation functions for
the polar, $\theta$ (middle row), and azimuthal, $\varphi$ (outermost right),
angles of vector $g_{2}$ $\rightarrow$P from POPC (top) and POPC/cholesterol
(bottom) CHARMM36 MD simulations. The grey lines represent data for each
individual phospholipid molecule in the simulation. $\theta$ and $\varphi$ are
the spherical angles in the laboratory coordinate frame defined by the
simulation box axes.  Figure S13: Orientation distribution, $p(\theta)$
(outermost left), and auto-correlation functions for the polar, $\theta$
(middle row), and azimuthal, $\varphi$ (outermost right), angles of vector
C1(sn-1)$\rightarrow$C1(sn-2) from POPC (top) and POPC/cholesterol (bottom)
CHARMM36 MD simulations. The grey lines represent data for each individual
phospholipid molecule in the simulation. $\theta$ and $\varphi$ are the
spherical angles in the laboratory coordinate frame defined by the simulation
box axes.


11institutetext: Anatoli Ivanov 22institutetext: Pennsylvania State University
at Wilkes-Barre, 44 University Drive, Dallas, PA 18612, USA, 22email:
<EMAIL_ADDRESS>33institutetext: Sergiy Shelyag 44institutetext: Flinders
University, Tonsley Innovation District, 1284 South Rd, Tonsley, SA, 5042,
Australia, 44email<EMAIL_ADDRESS>
# Explicit Periodic Solutions in a Delay Differential Equation

Anatoli Ivanov and Sergiy Shelyag

###### Abstract

We construct stable periodic solutions for a simple form nonlinear delay
differential equation (DDE) with a periodic coefficient. The equation involves
one underlying nonlinearity with the multiplicative periodic coefficient. The
well-known idea of reduction to interval maps is used in the case under
consideration, when both the defining nonlinearity and the periodic
coefficient are piece-wise constant functions. The stable periodic dynamics
persist under a smoothing procedure in a small neighborhood of the
discontinuity set. This work continues the research in recent paper IvaShe23
on stable periodic solutions of differential delay equations with periodic
coefficients.

##  1 Introduction

Scalar delay differential equations of the form

$x^{\prime}(t)=-\mu x(t)+f(x(t-\tau))$ (1)

are used as mathematical models in a broad variety of applications. Some of
the most well-known ones are the Mackey-Glass physiological models MacGla77 ,
Nicholson’s blowflies models GurBlyNis80 ; PerMalCou78 , and many others
GlaMac88 ; Kua93 . In spite of its formal simplicity equation (1) can exhibit
quite complex dynamical behaviors. Unlike scalar ordinary differential
equations the DDE (1) can have complicated dynamics of solutions including
oscillations about equilibrium, stable and unstable periodicity, and
complex/chaotic behaviors DieSvGSVLWal95 ; HalSVL93 .

More accurate mathematical models of real world phenomena take into account
some periodicity factors such as seasonal changes or circadian rhythm
parameter fluctuations and several others. Such improved models with respect
to equation (1) lead to similar type DDEs with periodic coefficients:

$x^{\prime}(t)=-\mu(t)x(t)+a(t)f(x(t-\tau)),$ (2)

where $\mu(t)$ and $a(t)$ are periodic functions with the same period. One of
the natural questions that appears in this context is whether DDE (2) has
periodic solutions with the same period as the periodic coefficients. This
work is an attempt to answer such question in a particular case.

Note that the periodicity problem we consider is different from similar ones
in many other publications, such as in e.g. Far17  (also see additional
references therein). We assume the negative feedback conditions on $f,x\cdot
f(x)<0,\;\forall x\neq 0,$ implying that DDE (2) admits the trivial solutions
$x\equiv 0$ (for any functions $\mu$ and $a$). This situation is similar to
the case of autonomous equation (1), when it admits a constant solution (with
the negative feedback in many cases). In Far17  and elsewhere, when such
constant positive equilibrium is perturbed by a multiplicative periodic
coefficient, the new problem does not admit constant solutions (positive
equilibrium) any longer.

A particular case of equation (2) is considered when $\mu(t)\equiv 0$ (see
equation (3) below). The initial choice of the nonlinearity $f$ and the
periodic coefficient $a(t)$ is as piecewise constant functions. Such simple
form allows for an explicit calculation of solutions for all forward times as
piecewise affine functions (made up continuously of line segments). The
problem of existence of slowly oscillating solutions and their stability is
reduced to the existence of fixed points of a simple interval map and their
attractivity. The dynamics of both is shown to persist when the two
discontinuous functions are replaced by close to them continuous functions (or
even of $C^{\infty}$ class).

The results of this paper are a continuation and an expansion of recently
obtained similar results by the same authors IvaShe23 .

##  2 Preliminaries

As in IvaShe23  we consider the special case of DDE (2) when $\mu(t)\equiv 0$:

$x^{\prime}(t)=a(t)f(x(t-1))$ (3)

Here initially the nonlinear function $f$ is given by
$f(x)=f_{0}(x)=-\mathrm{sign}(x)$ and the periodic coefficient
$a=a(t,a_{1},a_{2},p_{1},p_{2})$ is defined as

$a(t)=a_{0}(t)=\begin{cases}a_{1},\;\text{if}\;0\leq t<p_{1}\\\
a_{2},\;\text{if}\;p_{1}\leq t<p_{1}+p_{2}\\\ \text{periodic extension
outside}\;[0,p_{1}+p_{2})\;\text{for all}\;t\in\mathbf{R},\end{cases}$ (4)

with all the values $a_{1},a_{2},p_{1},p_{2}$ being positive and
$p_{1}+p_{2}:=T>1$.

For arbitrary initial function $\varphi(s)\in C=C([-1,0],\mathbf{R})$ such
that $\varphi>0\;\forall s\in[-1,0]$ the corresponding unique solution
$x=x(t,\varphi)$ exists for all $t>0$ (it is calculated by the consecutive
integration, the ”step method”). The solution is a continuous piece-wise
affine function made up of straight line segments. It is differentiable
everywhere except a countable set of points where two line segments with
different slopes match.

##  3 Main Results

###  3.1 Explicit Periodic Solutions

Given $f=f_{0}$ and $a=a_{0}$ as above for arbitrary initial function $\phi\in
C$ the corresponding solution $x=x(t,\phi)$ is constructed for all forward
times $t\geq 0$. It is a continuous piecewise differentiable function composed
of consecutive segments of affine (linear) functions.

We assume that periodic solutions have the form as depicted in Fig. 1 (this
will be shown to be the case for particular values of the parameters
$a_{1},a_{2},p_{1},p_{2}$). For arbitrary initial function $\phi\in C$ such
that $\phi(s)>0\;\forall s\in[-1,0],$ the corresponding forward solution
$x(t),t\geq 0,$ depends only on the value $\phi(0):=h>0$ and does not depend
on the other values $\phi(s),s\in[-1,0)$ of the initial interval. There exist
consecutive zeros $t_{1}>0$ and $t_{2}=t_{1}+2$ of the solution such that
$x(t)>0,t\in[0,t_{1})$, $x(t)<0,t\in(t_{1},t_{1}+2)$, and
$x(t)>0,t\in[t_{1}+2,p_{1}]$. The necessary condition for such first two zeros
to exist is that $p_{1}>2$. We also assume that $p_{1}<t_{1}+3<p_{1}+p_{2}$
and that the solution remains positive on the interval $[p_{1},p_{1}+p_{2}]$.
For the latter to be true one only has to require that
$x_{3}:=x(p_{1}+p_{2})>0$.

Figure 1: Slowly oscillating piece-wise affine solution

We next calculate the explicit values of
$t_{1},x_{1}=x(p_{1}),x_{2}=x(t_{1}+3),x_{3}=x(p_{1}+p_{2})$ in terms of the
parameters $a_{1},a_{2},p_{1},p_{2}$ and given initial value $h>0$. It is
straightforward to find that:

$t_{1}=\frac{h}{a_{1}},\quad x_{1}=x(p_{1})=-h+a_{1}p_{1}-2a_{1},\;$
$x_{2}=x(t_{1}+3)=\left(\frac{a_{2}}{a_{1}}-1\right)h+a_{1}p_{1}-2a_{1}+3a_{2}-a_{2}p_{1},$
$x_{3}=x(p_{1}+p_{2})=\left(2\frac{a_{2}}{a_{1}}-1\right)h+a_{1}(p_{1}-2)+a_{2}[6-(2p_{1}+p_{2})]:=mh+b$

In addition to the previous assumptions the condition $b>0$, i.e.

$b=a_{1}(p_{1}-2)+a_{2}[6-(2p_{1}+p_{2})]>0,$ (5)

guaranties that $x_{3}>0$, and thus that the solution
$x(t),t\in[0,p_{1}+p_{2}],$ is of desired shape (as shown in Fig. 1).

A fixed point $h_{*}>0$ of the affine map $F(x)=mx+b$ gives rise to a slowly
oscillating periodic solution $x=x(t,h_{*})$ of equation (1). Moreover, by the
linearity of $F$ and the construction, such periodic solution is
asymptotically stable if $|m|<1$. The latter is equivalent to $0<a_{2}<a_{1}$.
The unique fixed point $x=h_{*}$ is easily found as

$h_{*}=b/(1-m)=\frac{a_{1}[a_{1}(p_{1}-2)+a_{2}[6-(2p_{1}+p_{2})]]}{2(a_{1}-a_{2})}.$
(6)

Therefore, we arrive at the following statement

###### Theorem 3.1

Suppose that the parameters $a_{1},a_{2},p_{1},p_{2}$ are such that the
inequality (5) is satisfied and $a_{1}>a_{2}$. Then DDE (3) has an
asymptotically stable slowly oscillating periodic solution. The periodic
solution is generated by the initial function $\phi(s)\equiv
h_{*},s\in[-1,0],$ where $h_{*}$ is given by (6).

Note that due to the symmetry property of the nonlinearity $f$ (oddness) and
the procedure of construction of the periodic solution $x_{h_{*}}$, under the
assumption of Theorem 3.1 there also exists the symmetric to $x_{h_{*}}$
periodic solution generated by the initial function $\phi(s)\equiv-h_{*}$. The
two periodic solutions are related by $x_{-h_{*}}(t)\equiv-x_{h_{*}}(t)$.

###  3.2 Smoothing of Nonlinearities

In this subsection we demonstrate that the stable periodic solution derived in
the previous subsection persists under the standard smoothing procedure for
functions $f$ and $a$. The procedure is a well known one which has been used
in many publications (see e.g. IvaSha91 ; Pet80  and further references
therein). It consists in replacing each jump discontinuity by a continuous
affine function in a small neighborhood of a discontinuity point. We repeat
here the exact constructions as in IvaShe23 .

Consider the functions $f$ and $a$ as defined in subsection 3.1 :
$f(x)=f_{0}(x)$ and $a(t)=a_{0}(t)$. Let $\delta_{0}>0$ be small, and for
every $\delta\in(0,\delta_{0}]$ introduce the continuous functions
$f_{\delta}(x)$ and $a_{0}^{\delta}(t)$ by:

$f(x)=f_{\delta}(x)=\begin{cases}+1\;\text{if}\;x\leq-\delta\\\
-1\;\text{if}\;x\geq\delta\\\
-({1}/{\delta})x\;\text{if}\;x\in[-\delta,\delta],\end{cases}$ (7)

and

$a(t)=a_{\delta}(t)=\begin{cases}a_{2}+\frac{a_{1}-a_{2}}{2\delta}(t+\delta)\;\text{if}\;t\in[-\delta,\delta]\\\
a_{1}\;\text{if}\;t\in[\delta,p_{1}-\delta)\\\
a_{1}+\frac{a_{2}-a_{1}}{2\delta}[t-(p_{1}-\delta)]\;\text{if}\;t\in[p_{1}-\delta,p_{1}+\delta]\\\
a_{2}\;\text{if}\;t\in[p_{1}+\delta,p_{1}+p_{2}-\delta)\\\
a_{2}+\frac{a_{1}-a_{2}}{2\delta}[t-(p_{2}-\delta)]\;\text{if}\;t\in[p_{1}+p_{2}-\delta,p_{1}+p_{2}+\delta]\\\
\text{periodic extension on}\;\mathbf{R}\;\text{outside
interval}\;[0,p_{1}+p_{2}).\end{cases}$ (8)

The following statement is an extension of the main Theorem 3.1 to the case of
continuous functions $f_{\delta}$ and $a_{\delta}$.

###### Theorem 3.2

Suppose the assumptions of Theorem 3.1 are satisfied. There exists
$\delta_{0}>0$ such that for every $\delta\in[0,\delta_{0}]$ DDE (3) with
$f=f_{\delta}$ and $a=a_{\delta}$ has an asymptotically stable slowly
oscillating periodic solution $x_{\delta}(t).$ The solution $x_{\delta}(t)$ is
close in the uniform metric to the periodic solution $x_{h_{*}}$ of Theorem
3.1 for small $\delta\geq 0$, moreover $\lim_{\delta\to
0}\\{\sup_{t\in\mathbf{R}}|x_{\delta}(t)-x_{h_{*}}(t)|\\}=0$.

The proof is based on the explicit calculation of the similar value $x_{3}$ of
the solution $x(t,h)$ of the corresponding DDE (3). Its shape is similar to
that shown in Fig. 1, except that the corner points (where the derivative
$x^{\prime}$ is not defined) are replaced by a smooth $C^{1}$ curve, based on
the affine representation of both $f$ and $a$ around their discontinuity
points. The exact calculations repeat those in IvaShe23 . The final result is
a representation of $x_{3}(h)$ in the form $x_{3}(h)=\tilde{F}(h)$, where
$\tilde{F}(h)$ is $C^{1}$-close to the $F$ defined in subsection 3.1 (see
pages 74-75 of IvaShe23  for more details). Therefore, the existence,
stability, and closeness follow.

###  3.3 Numerical Demonstration

Periodic solutions described in subsections 3.1 and 3.2 can be easily verified
numerically. It is clear that sets of parameters leading to the existence of a
stable periodic solution, as described by Theorem 3.1, are not exceptional.
Any sufficiently small perturbation of a particular quadruple of the
parameters will produce an asymptotically stable periodic solution. They are
open in the set of all parameters. The exemplary parametric values for which
numerical periodic solutions have been obtained and confirmed to be of the
described form are given in Table 1.

$~{}~{}~{}~{}~{}a_{1}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}~{}a_{2}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}~{}p_{1}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}~{}p_{2}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}~{}h_{*}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}~{}T~{}~{}~{}~{}~{}$
---|---|---|---|---|---
1 | 0.25 | 2.5 | 1.5 | 0.25 | 4
2 | 0.5 | 2.5 | 2 | 1/3 | 4.5
2 | 0.25 | 2.5 | 1 | 4/7 | 3.5
1 | 0.5 | 3 | 1 | 0.5 | 4
2 | 1 | 3 | 1.5 | 0.5 | 4.5
2.5 | 0.5 | 3 | 4 | 0.31 | 7
3 | 0.5 | 3 | 4.5 | 0.45 | 7.5
5 | 0.5 | 3 | 3 | 1.94 | 6
5 | 1 | 3 | 2 | 1.88 | 5
Table 1: Parametric values for which stable periodic solutions of equation (3)
with the coefficient defined in equation (4) have been demonstrated
numerically.

##  4 Discussion

This paper demonstrates the existence of another type of stable periodic
solutions in DDE (3). They are complimentary to those which existence was
shown in our recent work IvaShe23 . There is a possibility of existence of
additional types of stable periodic solutions in (3) with piece-wise constant
defining functions $f$ and $a$, the question we intend to tackle in our next
research considerations. It is of interest to show the existence of stable
periodic solutions in a general DDE of the form (3), where $f$ and $a$ are
smooth functions of arbitrary type. An interesting and seemingly very
nontrivial problem is the existence of (stable) periodic solutions in DDE (3)
with smooth functions, when $f$ is arbitrary and fixed nonlinearity with the
negative feedback property and the periodic coefficient $a$ is parameter
dependent, e.g., $a=a(t,\varepsilon)$, such that $a(t,0)=a_{0}>0$ and
$a(t+T,\varepsilon)=a(t,\varepsilon)$, for some $T>0$. As the value
$\varepsilon=0$ gives periodic solutions in the well-known autonomous case of
equation (1), and some larger positive values of $\varepsilon_{0}$ result in
$T$-periodic solutions with the same period as $a$, the intermediate values of
$\varepsilon\in(0,\varepsilon_{0})$ may yield varied dynamics and complex
transition patterns. This is an open problem worthy further investigation.

###### Acknowledgements.

The authors thank the mathematical research institute MATRIX in Australia
where part of this research was performed. Its final version resulted from
collaborative activities of the authors during the workshop ”Delay
Differential Equations and Their Applications” (https://www.matrix-
inst.org.au/events/delay-differential-equations-and-their-applications/) held
in December 2023. The authors are also grateful for the financial support
provided for these research activities by Simons Foundation (USA), Flinders
University (Australia), and the Pennsylvania State University (USA). A.I.’s
research was also supported in part by the Alexander von Humboldt Stiftung
(Germany) during his visit to Justus-Liebig-Universität, Giessen, in June-
August 2023.

## References

  * (1) Diekmann, O., van Gils, S., Verdyn Lunel, S., Walther, H.-O.: Delay Equations: Complex, Functional, and Nonlinear Analysis. Springer-Verlag, Ser.: Applied Mathematical Sciences, Vol. 110 (1995).
  * (2) Faria, T. (2017). Periodic solutions for a non-monotone family of delayed differential equations with applications to Nicholson systems. J. Differential Equations, 263, no. 1, 509–533.
  * (3) Glass, L., Mackey, M.C.: From Clocks to Chaos. The Rhythms of Life. Prinston University Press, (1988).
  * (4) Gurney, W.S.C., Blythe, S.P., Nisbet, R.M. (1980). Nicholson’s blowflies revisited. Nature, 287, 17–21.
  * (5) Hale, J.K., Verduyn Lunel, S.M.: Introduction to Functional Differential Equations. Springer Applied Mathematical Sciences, vol. 99 (1993).
  * (6) Ivanov, A.F., Sharkovsky, A.N. (1991). Oscillations in singularly perturbed delay equations. Dynamics Reported (New Series), Springer Verlag, vol. 1, 165–224.
  * (7) Ivanov, A., Shelyag, S. (2023). Stable periodic solutions in scalar periodic differential delay equations. Archivum Mathemticum (Brno), 59 (1), 69–76.
  * (8) Kuang, Y.: Delay Differential Equations with Applications in Population Dynamics. Academic Press Inc. Series: Mathematics in Science and Engineering, Vol. 191, 398 pp. (1993).
  * (9) Mackey, M.C., Glass, L. (1977). Oscillation and chaos in physiological control systems. Science, 197, 287–289.
  * (10) Perez, J.F., Malta, C.P., Coutinho, F.A.B. (1978). Qualitative analysis of oscillations in isolated populations of flies. J. Theoret. Biol., 71, 505–514.
  * (11) Peters, H. (1980). Comportement chaotique d’une equation differentielle retardee. C. R. Acad. Sci. Paris Sér. A-B, 290, no. 24, A1119–-A1122. (in French)