Department of Biochemistry, University of Oxford, South Parks Road, Oxford,
OX1 3QU, United Kingdom

# Network Topology in Water Nanoconfined between Phospholipid Membranes

Fausto Martelli  IBM Research Europe, Hartree Centre, Daresbury, WA4 4AD,
United Kingdom<EMAIL_ADDRESS>Jason Crain  IBM Research Europe,
Hartree Centre, Daresbury, WA4 4AD, United Kingdom  Giancarlo Franzese  Secció
de Física Estadística i Interdisciplinària–Departament de Física de la Matèria
Condensada, Universitat de Barcelona, & Institut de Nanociència i
Nanotecnologia (IN2UB), Universitat de Barcelona, C. Martí i Franquès 1, 08028
Barcelona, Spain

###### Abstract

Water provides the driving force for the assembly and stability of many
cellular components. Despite its impact on biological functions, a nanoscale
understanding of the relationship between its structure and dynamics under
soft confinement has remained elusive. As expected, water in contact with
biological membranes recovers its bulk density and dynamics at $\sim 1$ nm
from phospholipid headgroups but surprisingly enhances its intermediate-range
order (IRO) over a distance, at least, twice as large. Here, we explore how
the IRO is related to the water’s hydrogen bond network (HBN) and its
coordination defects. We characterize the increased IRO by an alteration of
the HBN up to more than eight coordination shells of hydration water. The HBN
analysis emphasizes the existence of a bound-unbound water interface at $\sim
0.8$ nm from the membrane. The unbound water has a distribution of defects
intermediate between bound and bulk water, but with density and dynamics
similar to bulk, while bound water has reduced thermal energy and much more
HBN defects than low-temperature water. This observation could be fundamental
for developing nanoscale models of biological interactions and for
understanding how alteration of the water structure and topology, for example,
due to changes in extracellular ions concentration, could affect diseases and
signaling. More generally, it gives us a different perspective to study
nanoconfined water.

###### keywords:

water confined, phospholipid membrane, hydrogen bond, hydrogen bond network,
coordination defects, order parameter

It has been long recognized that the structure and function of biological
membranes are largely determined by the properties of hydration water, i.e.,
of the water in contact with the membrane 1. Indeed, the presence of water
strongly influences membrane stability, fluidity, and phase behavior, thereby
affecting membrane function and properties. Also, hydration water mediates the
interactions of biological membranes with other biomolecules and with ions 2,
3.

Biological membranes are composed of a large number of components, including
proteins, cholesterol, glycolipids, and ion channels, among others, but their
framework is provided by phospholipid molecules that self-assemble into
bilayers driven by the hydrophobic effect 4, 2, 3, 5, 6, 7, 8, 9, 10, 11, 12,
1, 13, 14, 15, 16. To decrease the interfacial free-energy, polar head groups
form contacts with water, while the apolar hydrocarbon tails minimize exposure
to water forming extended bilayers. Water is abundant in the interfacial
region of bilayers (lipid headgroups), establishing strong hydrogen bonds
(HBs) with the membrane. As a result of these strong interactions, the
orientational and translational dynamics of interfacial water is markedly
slowed down 7, 8, 9, 14, 15, 16, 6. Such slowing down has been observed also
in water in contact with proteins and sugars 17, 18, 19, 20, 21.

Recently we found an increase in the structural order at the intermediate
range when the dynamics of water confined by phospholipid membranes slows
down. This intermediate range order (IRO) propagates as far as (at least)
$\sim 2.4$ nm from the membrane surface 22, a larger distance than previously
calculated using other observables such as density and dynamical properties.
We recovered water’s bulk density and dynamical properties at a distance of
$\sim 1.2$ nm from the membrane surface 22, 23. Nonetheless, water is a
complex network-forming material, with a directional HBs network (HBN), the
topology of which is correlated to anomalous behavior, as we showed recently
24. Therefore, understanding how the HBN of water is affected by the
interactions with the membrane is of primary importance in understanding
biological properties at the molecular scale also in conditions in which a
biological membrane interacts with alien components such as, e.g., viruses.

In this article, we investigate the properties of water confined by
phospholipid membranes. Specifically, we measure the extent to which
phospholipid membranes affect the structural properties of water as well as
its HBN. As a typical model membrane, we use 1,2-Dimyristoyl-sn-
glycero-3-phosphocholine (DMPC) lipids. The DMPC is a phospholipid with a
choline headgroup and a tailgroup formed of two myristoyl chains(see fig. 1).
Choline-based phospholipids are ubiquitous in cell membranes and commonly used
in drug-targeting liposomes 25.

Figure 1:  Chemical structure of the DMPC phospholipid.

In this investigation, we probe the structural properties of water by
inspecting the IRO using a sensitive local metric, recently introduced by
Martelli et al. 26 and already applied in a wide variety of studies 26, 24,
22, 23, 27, 28. Next, we examine the topology and the quality of the HBN and
explore correlations of these observables with the behaviour of the structural
order. We probe the topology of the HBN via ring statistics, a tool widely
adopted in network-forming materials as a measure of closed loops present in
the network. Ring statistics have been previously employed in water to
characterize its different phases 26, 27, 29, 30, 31, as well as to
investigate the origin of water anomalies 24, 32, 33, 34, 35.

## Results

We inspect the IRO by computing the score function (Eq. 2), which provides a
nonlinear measure of the deviation in atomic or molecular arrangements from
those of a reference structure that, usually, corresponds to the medium’s
ground state at $T=0$ K. Therefore, we here compute $S$ using, as a reference,
the position of the oxygen in the second coordination shell in cubic ice,
i.e., a cuboctahedron ($\bar{C}$) that belongs to the class of Archimedean
solids enriched with edge transitivity 36 (inset in fig. 2). Similar results
hold if we use, as a reference, the position of the oxygen in the second
coordination shell in hexagonal ice.

At fixed $T=303$ K and $P=1$ atm, we compare the distributions of
$S_{\bar{C}}$ for bulk water and confined water at a distance of $1.6$ nm
$\leq z\leq 2.4$ nm (bin centered at $2.0$ nm) and $2.4$ nm $\leq z\leq 3.2$
nm (bin centered at $2.8$ nm) from the membrane (fig. 2 a). We emphasize that
the bin width of $0.8$ nm adopted in this work ensures that water molecules
centered in the middle of the bin have a second shell of neighbours falling
inside the same bin.

The distribution at $2.0$ nm does not match that of bulk water. In particular,
we find that water $2.0$ nm away from the membrane is more structured than
bulk water, with a $\sim 8\%$ increase of the $S_{\bar{C}}$ with maximum
probability, and a higher population in the large-$S_{\bar{C}}$ tail of the
probability distribution $P(S_{\bar{C}})$ (fig. 2 b).

On the other hand, the $P(S_{\bar{C}})$ for water at $2.8$ nm from the
membrane overlaps with that of bulk water, within our resolution, showing that
the value of the local order metric (LOM) of water between $2.4$ nm and $3.2$
nm is not affected by the membrane. Hence, our result implies that the effect
of the membrane on the structural properties of water should extend as far as
$1.6$ nm plus the second coordination shell distance ($\sim 0.45$ nm 37), and
at about $2.4$ nm minus $0.45$ nm, i.e., up to $(2.0\pm 0.05)$ nm,
approximately (fig. 2 b).

Since the properties of water emanate from the underlying network of HBs 24, a
natural question follows: Is there a connection between i) the observed
perturbations on the IRO of confined water, ii) the underlying HBN, and iii)
its quality in terms of broken and intact HBs? We will address this question
in the following discussion.

Figure 2:  Panel a): Probability distributions of score function $S_{\bar{C}}$
computed using the reference in the inset: the black line is for bulk water,
the green line is for water in the bin centered at $z=2.0$ nm, the red dashed
line is for water in the bin centered at $z=2.8$ nm. Inset: Reference made of
the oxygen positions (blue spheres) of the second coordination shell in cubic
ice. Blue sticks are a guide to the eyes to emphasize the geometrical
structure. Panel b): Difference $\Delta P$ between $P(S_{\bar{C}})$ for bulk
water and for the bin at $z=2.0$ nm (blue, continuous line), and between bulk
water and the bin at $z=2.8$ nm (orange, dashed line). Figure 3:  Probability
of the HB $n$-member rings, $P(n)$, computed from structures in bulk water
(open orange triangles), and in water in bins centered at different distances
from the membrane surfaces: $0.4$ nm (black dots), $1.2$ nm (red squares),
$2.0$ nm (green diamonds), and $2.8$ nm (blue triangles). All $P(n)$ are
normalized to unity and, therefore, do not reflect the total number of rings
of a given size.

We probe the HBN of water using the ring statistics, and we inspect the
quality of the network quantifying and characterizing coordination defects. We
compare the probability $P(n)$ of having an $n$-member ring,
$n\in\left[3,12\right]$, for the four bins that discretize the simulation box,
and for bulk water at the same thermodynamic conditions, $T=303$ K and $P=1$
atm (fig. 3).

For bulk water, as expected in diffusive media, $P(n)$ is broad and accounts
for very large rings. We find that for distances within the bin closer to the
bilayer, i.e., at a distance $z\leq 0.8$ nm from the membrane, $P(n)$ strongly
deviates from the corresponding probability in bulk water. Namely, we observe
a depletion in the number of larger rings and an increase, notably sharp in
the case of $n=6$, of shorter rings. We attribute the depletion of larger
rings to the proximity of the membrane, which represents a reduction of
dimensionality in the connectivity search pathways of water molecules. On the
other hand, we remark that $n=6$ represents the typical connectivity in
crystalline ice. Therefore, the increased number of hexagonal rings at $z=0.8$
nm indicates that, closer to the interface, the HBN seems to acquire a
topology closer to that of an ordered crystalline network. However, as we will
discuss later, this similarity is only apparent.

The increased number of short rings at $z\leq 0.8$ nm from the membrane is in
agreement with the dynamical slowing down 38, and the increment in the IRO 22
reported for similar distances. Hence, 1) the diffusion and rotational slowing
down 38, 2) the increased value of $S_{\bar{C}}$ (fig. 2), and 3) the
increased fraction of hexagonal rings (fig. 3) for water at $z\leq 0.8$ nm
from the membrane, suggest a connection between dynamics and structure as
measured by i) the positions of the oxygen atoms and, ii) the topology of the
HBN.

Moving from $z\leq 0.8$ nm to $0.8$ nm $<z\leq 2.4$ nm, we observe a marked
change in the distribution of $n$-rings with a decrease for $n\leqslant 6$ and
an increase for $n>6$ (fig. 3). We attribute the larger probability for
extended rings for $z>0.8$ nm to the increased dimensionality of the space
available. In particular, the difference in the $P(n)$ and in the
$S_{\bar{C}}$ (fig. 2), between bulk and water within the bin centered at
$z=2.0$ nm from the membrane, further points toward a close correlation
between structural properties at the level of the medium range, and the
topology of the HBN.

The drastic change in the ring probability between the bin centered in $z=0.4$
nm and the bins at a larger distance, is consistent with the recent discovery
of an interface between bound and unbound hydration water at about $0.5$ nm
from the membrane38. In Ref. 38 Calero and Franzese identify the interface
between i) the first hydration shell, partially made of water bound to the
membrane, with a structural role and an extremely slow dynamics, and ii) the
next shells with no water-lipids HBs and a dynamics ten time faster than bound
water, but still one order of magnitude slower than bulk water. Therefore,
ring probability can mark the structural difference between bound and unbound
water.

Moving to a distance $2.4$ nm $<z\leq 3.2$ nm from the surface, the $P(n)$
overlaps perfectly with the bulk case. Moreover, the $P(S_{\bar{C}})$ computed
within this bin (fig. 2) overlaps with the $P(S_{\bar{C}})$ of bulk water.
Therefore, we conclude that water recovers the structural (both IRO and HBN)
properties of bulk water only if at a distance larger than $2.4$ nm from the
membrane. This value is twice the $1.2$ nm at which water retrieves bulk
density and dynamics 22, 23.

To get further insights into the network topology, we inspect its quality.
When water is in the glass state, we can map its HBN to a nearly-hyperuniform
network, i.e., to a continuous random network characterized by a low fraction
of coordination defects and a suppression of long-range density fluctuations
39. Therefore, the number of broken HBs is a measure of the quality of the
HBN. In particular, it quantifies how far the HBN is from the two extreme
cases: a) the liquid and b) the continuous random network. Furthermore,
coordination defects directly affect the fluidity of liquid water. Therefore,
they can be related to water dynamics 40.

We perform a decomposition of the HBs per water molecule into acceptor-(A) and
donor-(D) types. We label as $\textit{A}_{2}\textit{D}_{2}$ a water molecule
with perfect coordination, i.e., donating two bonds and accepting two bonds.
We evaluate the quality of the HBN by computing the ratio of water molecules
that have different coordination, i.e., are not in the
$\textit{A}_{2}\textit{D}_{2}$ configuration. In particular, we focus our
attention on the following coordination configurations:
$\textit{A}_{1}\textit{D}_{1}$, $\textit{A}_{2}\textit{D}_{1}$,
$\textit{A}_{1}\textit{D}_{2}$, $\textit{A}_{2}\textit{D}_{2}$ and
$\textit{A}_{3}\textit{D}_{2}$. Other configurations do not contribute
significantly 41.

Figure 4:  Percentage-wise decomposition of the intact HBs per water molecule
into acceptor-(A) and donor-(D) for water in bins centered at different
distances from the membrane and for bulk water. Sets are for bins at $0.4$ nm
(black dots), $1.2$ nm (red squares), $2.0$ nm (green diamonds), $2.8$ nm
(blue triangles) and bulk (orange open triangles). The $x$-axis labels
$\textit{A}_{x}\textit{D}_{y}$ indicate the number of acceptor
($\textit{A}_{x}$) and donor ($\textit{D}_{y}$) HBs, respectively, of the
configurations schematically represented on the panel’s top (with the oxygen
of central water molecule in blue). For clarity we omit combinations with
minor contributions, e.g., $\textit{A}_{3}\textit{D}_{1}$,
$\textit{A}_{0}\textit{D}_{y}$, $\textit{A}_{x}\textit{D}_{0}$, etc.

We compare the percentage of intact HBs for bulk and confined water, as a
function of the distance from the membrane (fig. 4). We find that the HBN in
bulk water is dominated by $\textit{A}_{2}\textit{D}_{2}$ ($\sim 37\%$)
perfect coordinations. Water molecules involved in three HBs in the form of
the defect $\textit{A}_{1}\textit{D}_{2}$ comprise the next largest percentage
($\sim 20\%$), followed by the $\textit{A}_{2}\textit{D}_{1}$ and
$\textit{A}_{1}\textit{D}_{1}$ types and, finally, by the
$\textit{A}_{3}\textit{D}_{2}$.

This result, based on TIP3P water, is in agreement with the trend in _ab
initio_ liquid water at ambient conditions examined with different functionals
41. In particular, in _ab initio_ liquid water, the frequency of
$\textit{A}_{1}\textit{D}_{2}$ is almost twice that of
$\textit{A}_{2}\textit{D}_{1}$ at all levels of theories 41.

Close to the surface of the membrane, at $z\leq 0.8$ nm, the network of HBs
largely deviates from that of bulk water. The network is dominated by
$\textit{A}_{1}\textit{D}_{1}$ and $\textit{A}_{1}\textit{D}_{2}$ defects
($\sim 25\%$), followed by $\textit{A}_{2}\textit{D}_{1}$ and
$\textit{A}_{2}\textit{D}_{2}$ configurations ($\sim 15\%$), and a small
percentage of higher coordination defects $\textit{A}_{3}\textit{D}_{2}$
($\sim 3\%$). Such composition is very consistent with the results found for
bound water at $z\leq 0.5$ nm 38. In particular, we find here the same
percentage of defects with three water-water HBs ($40\%$). Furthermore, we
observe numbers, very close to those in Ref.38 , for perfectly coordinated
configurations ($\sim 20\%$), and defects with two ($\sim 30\%$) and five
water-water HBs ($\sim 1\%$).

However, close to the membrane, the decrease of perfectly coordinated water
molecules seems to be inconsistent with the higher local order of water 22,
23, and also with the enhanced contribution of six-fold rings (fig. 3). This
discrepancy is only apparent, for two reasons. First, both the IRO and the
ring statistics are a measure of local order beyond short range, while the
quality of the HBN is strictly a short-range measure. Second, our calculations
include only water-water HBs and do not account for the (strong) HBs between
water molecules and the phospholipid headgroups 42, 43, 38. Instead, $\sim
30\%$ of the water molecules in the first hydration shell are bound to the
membrane with at least one HB 38. This observation explains why the dynamical
slowing down 22 of bound water can be interpreted as a local reduction of
thermal noise that allows water molecules to organize in space in more ordered
geometrical configurations 38.

Moving away from the surface, at a distance of $0.8$ nm$<z\leq 2.4$ nm, the
most appreciable effect on the quality of the HBN is a marked reduction of
$\textit{A}_{1}\textit{D}_{1}$ defects down to $\sim 18\%$, mostly accounting
for the absence of HBs between water molecules and phospholipid headgroups 38,
and a corresponding drastic increase in the percentage of perfectly
coordinated water molecules ($\textit{A}_{2}\textit{D}_{2}$) up to $\sim
25\%$, confirming the analysis done for unbound water 38.

At these distances bulk density and dynamical properties of water are
recovered almost fully 22, 38. However, the quality (defects) of the HBN
strongly deviates from that of bulk water, accounting for its different
topology (ring probability, fig. 3) and its different structural properties
22.

Upon increasing the distance from the membrane, we find a reduction of most of
the coordination defects and a corresponding increase of perfectly coordinated
water molecules, i.e., an improvement in the quality of the HBN (fig. 4).
Nevertheless, we recover the bulk-like composition only at a distance $z>2.4$
nm from the membrane, as in our analysis for both the IRO (fig. 2) and the
topology of the HBN (fig. 3).

It is interesting to observe that the percentage of the defect type
$\textit{A}_{2}\textit{D}_{1}$ is mostly constant in all bins and, therefore,
is independent of the distance from the membrane (fig. 4). We are currently
working on rationalizing this intriguing evidence.

## Conclusions

The relation between water dynamics and structure is elusive in bulk 44 and
even more under nanoconfinement 45, especially when the confining surfaces are
soft 46. On the other hand, the relationships among the hydration structure
and molecular fluidity at membrane/water interfaces are relevant in many
biological processes 47. Here, we study why water recovers its density and
dynamics at $\sim 1.2$ nm from a membrane 22, 38 while has an intermediate
range order (IRO) 26 higher than bulk up to a distance twice as large 22. To
understand this surprising result, we focus on the hydrogen bond network
(HBN), analyzing its topology (ring statistics) and its quality (population of
perfectly coordinated water molecules and defects). We find that the increased
IRO is characterized by an alteration of the HBN.

In particular, for bound water 38, i.e., water at short distances (here less
than $0.8$ nm) from the membrane, we show that the HBN topology and quality
are very different from those of low-temperature bulk water. Although bound
water has an HBN with a large fraction of hexagonal rings as in crystalline
water, it has a much higher number of defects than low-temperature water. We
find that $\textit{A}_{1}\textit{D}_{1}$ and $\textit{A}_{1}\textit{D}_{2}$
account together for 50% of all the defects due to water strong HBs with the
membrane. These strong HBs locally reduce the water’s thermal energy and slow
down its dynamics.

We show that the HBN analysis is able to mark the existence of the bound-
unbound water interface 38. We find a sudden qualitative change in the ring
statistics for hydration water at a distance $z>0.8$ nm from the membrane.
Also the defects distribution clearly shows that water in the range $0.8$
nm$<z\leq 2.4$ nm is neither bound to the membrane, neither bulk. Indeed, it
has much less $\textit{A}_{1}\textit{D}_{1}$ defects than bound water, and
much less perfectly-coordinated molecules than bulk. Nevertheless, at these
distances, the structural differences between unbound and bulk water are
disguised in water’s density and dynamics 38.

The difference in topology and defects smear out at distances larger than
$2.4$ nm. This distance corresponds to more than eight coordination shells of
hydration water. Hence, our results support the evidence of long-range effects
measured in terahertz and dielectric relaxation experiments 48, 9, 49, 50. We
expect our conclusions to hold and eventually be emphasized by water
potentials more realistic than TIP3P, which is quite poor in terms of
structural properties beyond the short-range.

Our findings should be taken into account when interpreting experimental
results and when developing membrane-water interaction potentials. They can
help in better understanding water in biological processes at large, in
particular those where hydration or structural changes play a role. Variations
of ions concentration drastically change the water HBN 51 and its dynamics 52,
with an effect that is similar to an increase of pressure 53, or a decrease of
temperature for dehydration 38. These variations in the extracellular matrix
can promote, for example, cardiac disease and arterial hardening in healthy
men 54 or atherosclerosis and inflammatory signaling in endothelial cells 55.
Hence, our results entail further investigation about the relationship between
this category of diseases with the water HBN rearrangements due to changes in
hydration or ionic concentrations.

## Methods

### Simulation details

The systems considered here have the same geometry as in our previous
simulations 22 but with a 15% increase in hydration, i.e., they are composed
of 128 DMPC lipids in a bilayer and $8100$ water molecules, with periodic
boundary conditions in such a way that water is confined between the two sides
of two replicas of the same membrane. We perform molecular dynamics (MD)
simulations on IBM POWER8 machines with NVIDIA Kepler K80 GPUs using the
simulation package NAMD 2.9 56 at a temperature of $T=303$ K and an average
pressure of $p=1$ atm. We set the simulation timestep to $2$ fs. We describe
the structure of phospholipids and their mutual interactions by the recently
parameterized force field CHARMM36 57, 58, which is able to reproduce the area
per lipid in excellent agreement with experimental data. The water model
employed in our simulations, consistent with the parametrization of CHARMM36,
is the modified TIP3P 59. We cut off the Van der Waals interactions at $12$ Å
with a smooth switching function starting at $10$ Å. We compute the long-
ranged electrostatic forces with the particle-mesh Ewald method 60, using a
grid spacing of $1$ Å. Our simulation box is anisotropic, with $L_{z}>L_{x}$,
with $L_{x}=L_{y}$. This anisotropy ensures that there are no errors caused by
the calculation of long-range electrostatics. During the $NpT$ simulations we
always keep this condition, with $\sim 5\%$ fluctuations for the values of
$L_{x}$, $L_{y}$, and $L_{z}$. After energy minimization, we equilibrate the
hydrated phospholipid bilayers for $10$ ns followed by a production run of $2$
ns in the $NpT$ ensemble at $p=1$ atm. The energy profile is shown in fig. 5.

Figure 5:  Energy profile for the system under consideration. The red dashed
arrow defines the limit after which we start the production and data analysis.

In the simulations, we control the temperature with a Langevin thermostat 61
using a damping coefficient of $0.1$ ps-1 and we control the pressure by a
Nosé-Hoover Langevin barostat 62 with a piston oscillation time of $0.2$ ps
and a damping time of $0.1$ ps. We also perform numerical simulations of bulk
TIP3P water (4000 molecules) adopting the same protocol and at the same
thermodynamic conditions. It is worthy to mention, at this point, that the
isotropic $NpT$ ensures that experimental observables such as, e.g., the area
per lipid and NMR order parameters, are properly reproduced 63, 64, 65.

In order to investigate the IRO and the HBN, we divide the systems along the
direction perpendicular to the phospholipid bilayer in equally spaced bins
such that the thickness of each bin is $0.8$ nm. This thickness is chosen in
such a way that a molecule of water in the bin’s center has, approximately,
it’s second coordination shell included in the same bin at $T=303$ K and $P=1$
atm. For our hydration and bin’s size, the two membranes are separated by
eight bins, hence we can analyze four different distances between 0 and $3.2$
nm. All the distances are measured taking as reference distance, for each side
of the membrane, the position where the phospholipid density distribution, at
thermodynamic equilibrium, has a maximum. We measure the observables of
interest for each water molecule within each bin centered at a distance $z$
from the center of the bilayer. It is worthy to mention that several more
sophisticated ways of computing the distances from a rough membrane surface
have been reported in the literature 38, 66. On the other hand, such methods
show differences whit respect to our approach when thinner bins are
implemented and only in the proximity of the surface.

### The local order metric

We here briefly discuss the basic ideas behind the LOM we have used to probe
the structural properties of water. Details can be found in Ref. 26.

The local environment of a water molecule $j$ in a snapshot defines a local
pattern formed by $M$ neighboring sites. Here, we consider only the oxygen
atoms’ second neighbors of the oxygen of the molecule $j$.

There are $N$ local patterns, one for each atomic site $j$ in the system.
Indicating by $\mathbf{P}_{i}^{j}(i=1,M)$ the position vectors in the
laboratory frame of the $M$ neighbors of site $j$, their centroid is given by
$\mathbf{P}_{c}^{j}\equiv\frac{1}{M}\sum_{i=1}^{M}\mathbf{P}_{i}^{j}$. In the
following we refer the positions of the sites of the pattern to their
centroid, i.e.
$\mathbf{P}_{i}^{j}-\mathbf{P}_{c}^{j}\rightarrow\mathbf{P}_{i}^{j}$.

The local reference is a set of $M$ sites, labeled by indices $i(i=1,M)$,
located at positions $\mathbf{R}_{i}^{j}$ around the molecules $j$ in ideal
positions, typically as in a lattice of choice. The step of the reference
lattice is fixed equal to equilibrium O–O distance, $d$, in the water
coordination shell at the thermodynamic conditions of interest.

For each oxygen site $j$ the centroid of the reference is set to coincide with
the centroid of the pattern. The reference orientation is, instead, arbitrary,
forming angles $\theta,\phi,\psi$ with the pattern.

The LOM $S(j)$ at site $j$ is the maximum of the overlap function with respect
to the orientation of the reference and the permutation of the pattern
indices,

$S(j)\equiv\max_{\theta,\phi,\psi;\mathcal{P}}\prod_{i=1}^{M}\exp\left(-\frac{\left|\mathbf{P}_{i_{\mathcal{P}}}^{j}-\mathbf{R}_{i}^{j}\right|^{2}}{2\sigma^{2}M}\right).$
(1)

Here $i_{\mathcal{P}}$ are the permuted indices of the pattern sites
corresponding to a permutation $\mathcal{P}$, and $\sigma=d/4.4$ is a
parameter that controls the spread of the Gaussian functions.

If $L$ is the number of proper point symmetry operations of the reference, the
overlap function (Eq. 1) has $L$ equivalent maxima. Therefore, it is
sufficient to compute $S(j)$ for only a fraction $1/L$ of the Euler angle
domain $\Omega$, which we may call $\Omega/L$, the irreducible domain of the
Euler angles. Inside $\Omega/L$ we pick at random, with uniform probability,
$15$ orientations and we optimize them using a conjugate gradients procedure.

The LOM is an intrinsic property of the local environment at variance with the
overlap function $\mathcal{O}(j)$ that depends on the orientation of the
reference and on the ordering of the sites in the pattern. The LOM satisfies
the inequalities $0\lesssim S(j)\leq 1$. The two limits correspond,
respectively, to a completely disordered local pattern ($S(j)\rightarrow 0$)
and to an ordered local pattern matching perfectly the reference
($S(j)\rightarrow 1$), therefore grading each local environment on an
increasing scale of local order from zero to one.

The order parameter score function $S$ is the site-averaged LOM:

$S\equiv\frac{1}{N}\sum_{j=1}^{N}S(j),$ (2)

### Definition of rings

Several definitions of rings and counting schemes have been reported in the
literature 67, 68, 69, 70, 71, 72, 73. Recently, Formanek and Martelli have
shown that different schemes allow us to access different information 74.

Here, we construct rings as in fig. 6. We adopt the geometric definition of HB
75, that is in qualitative agreement with other definitions over a wide range
of thermodynamic conditions 76, 77. We start from a tagged water molecule and
recursively traverse the HBN until we reached again the starting point, or we
exceed the maximal ring size considered, 12 water molecules in our case. We
consider only the primitive rings, i.e., rings that can not be decomposed into
smaller ones 78, 71, 67. As shown in Ref. 74, this definition provides a rich
amount of information about the network.

Figure 6:  Schematic representation of the ring construction for a given HBN
between water molecules. We start from water molecule labeled as 1 (O atoms
are represented in red, H atoms in white). By following the directional HBs
from H to O (blue arrows), we cross the HBN until we return to molecule 1 or
we exceeds 12 steps and then take only those rings that cannot be decomposed
in sub-rings. Here, we find a hexagonal ring from molecule 1 to 6 and and a
pentagonal ring from 1 to 9.

## Acknowledgements

F.M. and J.C. acknowledge support from the STFC Hartree Centre’s Innovation
Return on Research programme, funded by the Department for Business, Energy
and Industrial Strategy. G.F. acknowledges support from the Spanish grant
PGC2018-099277-B-C22 (MCIU/ AEI/ ERDF) and ICREA Foundation (ICREA Academia
prize).

## References

  * Berkowitz et al. 2006 Berkowitz, M. L.; Bostick, D. L.; Pandit, S. Aqueous Solutions next to Phospholipid Membrane Surfaces: Insights from Simulations. _Chem. Rev._ 2006, _106_ , 1527–1539
  * Yeagle 2011 Yeagle, P. L., Ed. _The Structure of Biological Membranes_ , 3rd ed.; CRC Press: New York, 2011
  * Berkowitz 2019 Berkowitz, M. L., Ed. _Biomembrane Simulations: Computational Studies of Biological Membranes_ , 1st ed.; Series in Computational Biophysics; CRC Press: New York, 2019
  * Fitter et al. 1999 Fitter, J.; Lechner, R. E.; Dencher, N. A. Interactions of Hydration Water and Biological Membranes Studied by Neutron Scattering. _J. Phys. Chem. B_ 1999, _103_ , 8036–8050
  * Trapp et al. 2010 Trapp, M.; Gutberlet, T.; Juranyi, F.; Unruh, T.; Demé, B.; Tehei, M.; Peters, J. Hydration Dependent Studies of Highly Aligned Multilayer Lipid Membranes by Neutron Scattering. _J. Chem. Phys._ 2010, _133_ , 164505
  * Wassall 1996 Wassall, S. R. Pulsed Field Gradient-Spin Echo NMR Studies of Water Diffusion in a Phospholipid Model Membrane. _Biophys. J._ 1996, _71_ , 2724–2732
  * Volkov et al. 2007 Volkov, V. V.; Palmer, D. J.; Righini, R. Distinct Water Species Confined at the Interface of a Phospholipid Membrane. _Phys. Rev. Lett._ 2007, _99_ , 078302
  * Zhao et al. 2008 Zhao, W.; Moilanen, D. E.; Fenn, E. E.; Fayer, M. D. Water at the Surfaces of Aligned Phospholipid Multibilayer Model Membranes Probed with Ultrafast Vibrational Spectroscopy. _J. Am. Chem. Soc._ 2008, _130_ , 13927–13937
  * Tielrooij et al. 2009 Tielrooij, K. J.; Paparo, D.; Piatkowski, L.; Bakker, H. J.; Bonn, M. Dielectric Relaxation Dynamics of Water in Model Membranes Probed by Terahertz Spectroscopy. _Biophys. J._ 2009, _97_ , 2848–2492
  * Hua et al. 2015 Hua, W.; Verreault, D.; Allen, H. C. Solvation of Calcium–Phosphate Headgroup Complexes at the DPPC/Aqueous Interface. _ChemPhysChem_ 2015, _16_ , 3910–3915
  * Róg et al. 2002 Róg, T.; Murzyn, K.; Pasenkiewicz-Gierula, M. The Dynamics of Water at the Phospholipid Bilayer Surface: A Molecular Dynamics Simulation Study. _Chem. Phys. Lett._ 2002, _352_ , 323–327
  * Bhide and Berkowitz 2005 Bhide, S. Y.; Berkowitz, M. L. Structure and Dynamics of Water at the Interface with Phospholipid Bilayers. _J. Chem. Phys._ 2005, _123_ , 224702
  * von Hansen et al. 2013 von Hansen, Y.; Gekle, S.; Netz, R. R. Anomalous Anisotropic Diffusion Dynamics of Hydration Water at Lipid Nembranes. _Phys. Rev. Lett._ 2013, _111_ , 118103
  * Zhang and Berkowitz 2009 Zhang, Z.; Berkowitz, M. L. Orientational Dynamics of Water in Phospholipid Bilayers with Different Hydration Levels. _J. Phys. Chem. B_ 2009, _113_ , 7676–7680
  * Gruenbaum and Skinner 2011 Gruenbaum, S. M.; Skinner, J. L. Vibrational Spectroscopy of Water in Hydrated Lipid Multi–Bilayers. I. Infrared Spectra and Ultrafast Pump–Probe Observables. _J. Chem. Phys._ 2011, _135_ , 075101
  * Borallo et al. 2016 Borallo, C. C.; Stanley, E. H.; Franzese, G. Structural Interpretation of the Large Slowdown of Water Dynamics at Stacked Phospholipid Membranes for Decreasing Hydration Level: All-Atom Molecular Dynamics. _Materials_ 2016, _9_ , 319
  * Camisasca et al. 2018 Camisasca, G.; Iorio, A.; Marzio, M. D.; Gallo, P. Structure and Slow Dynamics of Protein Hydration Water. _J. Mol. Liq._ 2018, _268_ , 903–910
  * Iorio et al. 2019 Iorio, A.; Camisasca, G.; Rovere, M.; Gallo, P. Characterization of Hydration Water in Supercooled Water–Trehalose Solutions: The Role of the Hydrogen Bonds Network. _J. Chem. Phys._ 2019, _51_ , 044507
  * Iorio et al. 2019 Iorio, A.; Camisasca, G.; Gallo, P. Glassy Dynamics of Water at Interface with Biomolecules: A Mode Coupling Theory Test. _Sci. China Phys. Mech._ 2019, _62_ , 107011
  * Iorio et al. 2019 Iorio, A.; Camisasca, G.; Gallo, P. Slow Dynamics of Hydration Water and the Trehalose Dynamical Transition. _J. Mol. Liq._ 2019, _282_ , 617–625
  * Iorio et al. 2020 Iorio, A.; Minozzi, M.; Camisasca, G.; Rovere, M.; Gallo, P. Slow Dynamics of Supercooled Hydration Water in Contact with Lysozyme: Examining the Cage Effect at Different Length Scales. _Philos. Mag_ 2020, _0_ , 1–14
  * Martelli et al. 2018 Martelli, F.; Ko, H.-Y.; Borallo, C. C.; Franzese, G. Structural Properties of Water Confined by Phospholipid Membranes. _Front. Phys._ 2018, _13_ , 136801
  * Samatas et al. 2019 Samatas, S.; Calero, C.; Martelli, F.; Franzese, G. In _Biomembrane Simulations Computational Studies of Biological Membranes_ ; Berkowitz, M., Ed.; CRC Press: New York, 2019; p 69
  * Martelli 2019 Martelli, F. Unravelling the Contribution of Local Structures to the Anomalies of Water: The Synergistic Action of Several Factors. _J. Chem. Phys._ 2019, _150_ , 094506
  * Hamley 2007 Hamley, W. _Introduction to Soft Matter_ ; John Wiley and Sons, West Sussex, England, 2007
  * Martelli et al. 2016 Martelli, F.; Ko, H.-Y.; Oğuz, E. C.; Car, R. Local-Order Metric for Condensed Phase Wnvironments. _Phys. Rev. B_ 2016, _97_ , 064105
  * Martelli et al. 2018 Martelli, F.; Giovambattista, N.; Torquato, S.; Car, R. Searching for Crystal-Ice Domains in Amorphous Ices. _Phys. Rev. Materials_ 2018, _2_ , 075601
  * Santra et al. 2018 Santra, B.; Ko, H.-Y.; Yeh, Y.-W.; Martelli, F.; Kaganovich, I.; Raitses, Y.; Car, R. Root-Growth of Boron Nitride Nanotubes: Experiments and Ab Initio Simulations. _Nanoscale_ 2018, _10_ , 22223
  * Marton̆ák et al. 2004 Marton̆ák, R.; Donadio, D.; Parrinello, M. Polyamorphism of Ice at Low Temperatures from Constant-Pressure Simulations. _Phys. Rev. Lett._ 2004, _92_ , 225702
  * Marton̆ák et al. 2005 Marton̆ák, R.; Donadio, D.; Parrinello, M. Evolution of the Structure of Amorphous Ice: From Low-Density Amorphous through High-Density Amorphous to Very High-Density Amorphous Ice. _J. Chem. Phys._ 2005, _122_ , 134501
  * Camisasca et al. 2019 Camisasca, G.; Schlesinger, D.; Zhovtobriukh, I.; Pitsevich, G.; Pettersson, L. G. M. A Proposal for the Structure of High– and Low–Density Fluctuations in Liquid Water. _J. Chem. Phys._ 2019, _151_ , 034508
  * Russo and Tanaka 2014 Russo, J.; Tanaka, H. Understanding Water’s Anomalies with Locally Favoured Structures. _Nat. Commun._ 2014, _5_ , 3556
  * Palmer et al. 2014 Palmer, J. C.; Martelli, F.; Liu, Y.; Car, R.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Metastable Liquid–Liquid Transition in a Molecular Model of Water. _Nature_ 2014, _510_ , 385–388
  * Santra et al. 2015 Santra, B.; DiStasio Jr., R. A.; Martelli, F.; Car, R. Local Structure Analysis in $AbInitio$ Liquid Water. _Mol. Phys._ 2015, _113_ , 2829–2841
  * Ansari et al. 2020 Ansari, N.; Ansari, N.; Onat, B.; Sosso, G. C.; Hassanali, A. Insights into the Emerging Networks of Voids in Simulated Supercooled Water. _J. Phys. Chem. B_ 2020, _124_ , 2180–2190
  * Read and Wilson 2016 Read, R. C.; Wilson, J. R. _An Atlas of Graphs_ ; Oxford University Press: Oxford, 2016
  * Mark and Nilsson 2001 Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. _J. Phys. Chem. A_ 2001, _105_ , 9954–9960
  * Calero and Franzese 2019 Calero, C.; Franzese, G. Membranes with Different Hydration Levels: The Interface between Bound and Unbound Hydration Water. _J. Mol. Liq._ 2019, _273_ , 488–496
  * Martelli et al. 2017 Martelli, F.; Torquato, S.; Giovambattista, N.; Car, R. Large-Scale Structure and Hyperuniformity of Amorphous Ices. _Phys. Rev. Lett._ 2017, _119_ , 136002
  * de los Santos and Franzese 2012 de los Santos, F.; Franzese, G. Relations between the Diffusion Anomaly and Cooperative Rearranging Regions in a Hydrophobically Nanoconfined Water Monolayer. _Phys. Rev. E_ 2012, _85_ , 010602–
  * DiStasio Jr. et al. 2014 DiStasio Jr., R. A.; Santra, B.; Li, Z.; Wu, X.; Car, R. The Individual and Collective Effects of Exact Exchange and Dispersion Interactions on the Ab Initio Structure of Liquid Water. _J. Chem. Phys._ 2014, _141_ , 084502
  * Binder 2003 Binder, H. The Molecular Architecture of Lipid Membranes: New Insights from Hydration–Huning Infrared Linear Dichroism Spectroscopy. _Appl. Spectrosc. Rev._ 2003, _30_ , 15–69
  * Chen et al. 2010 Chen, X.; Hua, W.; Huang, Z.; Allen, H. C. Interfacial Water Structure Associated with Phospholipid Membranes Studied by Phase–Sensitive Vibrational Sum Frequency Generation Spectroscopy. _J. Am. Chem. Soc._ 2010, _132_ , 1336–11342
  * Verde et al. 2019 Verde, A. R.; Montes de Oca, J. M.; Accordino, S. R.; Alarcón, L. M.; Appignanesi, G. A. Comparing the Performance of Two Structural Indicators for Different Water Models while Seeking for Connections between Structure and Dynamics in the Glassy Regime. _J. Chem. Phys._ 2019, _150_ , 244504
  * Joseph and Aluru 2008 Joseph, S.; Aluru, N. R. Why Are Carbon Nanotubes Fast Transporters of Water? _Nano Letters_ 2008, _8_ , 452–458
  * Ruiz Pestana et al. 2018 Ruiz Pestana, L.; Felberg, L. E.; Head-Gordon, T. Coexistence of Multilayered Phases of Confined Water: The Importance of Flexible Confining Surfaces. _ACS Nano_ 2018, _12_ , 448–454
  * Asakawa et al. 2012 Asakawa, H.; Yoshioka, S.; Nishimura, K.-i.; Fukuma, T. Spatial Distribution of Lipid Headgroups and Water Molecules at Membrane/Water Interfaces Visualized by Three-Dimensional Scanning Force Microscopy. _ACS Nano_ 2012, _6_ , 9013–9020
  * Ebbinghaus et al. 2007 Ebbinghaus, S.; Kim, S. J.; Heyden, M.; Yu, X.; Heugen, U.; Gruebele, M.; Leitner, D. M.; Havenith, M. An Extended Dynamical Hydration Shell around Proteins. _Proc. Natl. Ac. Sci. USA_ 2007, _104_ , 20749–20752
  * Zhang et al. 2013 Zhang, C.; Gygi, F.; Galli, G. Strongly Anisotropic Dielectric Relaxation of Water at the Nanoscale. _J. Phys. Chem. Lett._ 2013, _4_ , 2477–2481
  * Hishida and Tanaka 2013 Hishida, M.; Tanaka, K. Long–Range Hydration Effect of Lipid Membrane Studied by Terahertz Time-Domain Spectroscopy. _Phys. Rev. Lett._ 2013, _106_ , 158102
  * Mancinelli et al. 2007 Mancinelli, R.; Botti, A.; Bruni, F.; Ricci, M. A.; Soper, A. K. Perturbation of Water Structure Due to Monovalent Ions in Solution. _Phys. Chem. Chem. Phys._ 2007, _9_ , 2959–2967
  * Fayer et al. 2009 Fayer, M. D.; Moilanen, D. E.; Wong, D.; Rosenfeld, D. E.; Fenn, E. E.; Park, S. Water Dynamics in Salt Solutions Studied with Ultrafast Two-Dimensional Infrared (2D IR) Vibrational Echo Spectroscopy. _Acc. Chem. Res._ 2009, _42_ , 1210–9
  * Gallo et al. 2014 Gallo, P.; Corradini, D.; Rovere, M. Do Ions Affect the Structure of Water? The Case of Potassium Halides. _J. Mol. Liq._ 2014, _189_ , 52–56
  * Arnaoutis et al. 2017 Arnaoutis, G.; Kavouras, S. A.; Stratakis, N.; Likka, M.; Mitrakou, A.; Papamichael, C.; Sidossis, L. S.; Stamatelopoulos, K. The Effect of Hypohydration on Endothelial Function in Young Healthy Adults. _Eur. J. Nutr._ 2017, _56_ , 1211–1217
  * Dmitrieva and Burg 2015 Dmitrieva, N. I.; Burg, M. B. Elevated Sodium and Dehydration Stimulate Inflammatory Signaling in Endothelial Cells and Promote Atherosclerosis. _PLOS ONE_ 2015, _10_ , 1–22
  * Phillips et al. 2005 Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. _J. Comput. Chem._ 2005, _26_ , 1781–1802
  * Klauda et al. 2010 Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. _J. Phys. Chem. B_ 2010, _114_ , 7830–7843
  * Lim et al. 2012 Lim, J. B.; Rogaski, B.; Klauda, J. B. Update of the Cholesterol Force Field Parameters in CHARMM. _J. Phys. Chem. B_ 2012, _116_ , 203–210
  * Jorgensen et al. 1983 Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. Comparison of Simple Potential Functions for Simulating Liquid Water. _J. Chem. Phys._ 1983, _79_ , 926
  * Essmann et al. 1995 Essmann, U.; Perera, L.; Berkowitz, M. L. A Smooth Particle Mesh Ewald Method. _J. Chem. Phys._ 1995, _103_ , 8577
  * Berendsen et al. 1984 Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. _J. Phys. Chem._ 1984, _81_ , 3684
  * Feller et al. 1995 Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. _J. Phys. Chem._ 1995, _103_ , 4613
  * Taylor et al. 2009 Taylor, J.; Whiteford, N. E.; Bradley, G.; Watson, G. W. Validation of All–Atom Phosphatidylcholine Lipid Force Fields in the Tensionless NPT Ensemble. _Biochim. Biophys. Acta_ 2009, _1788_ , 638–649
  * Davis et al. 2009 Davis, J. E.; Rahaman, O.; Patel, S. Molecular Dynamics Simulations of a DMPC Bilayer Using Nonadditive Interaction Nodels. _Biophys. J.,_ 2009, _96_ , 285–402
  * Gapsys et al. 2013 Gapsys, V.; de Groot, B. L.; Briones, R. Computational Analysis of Local Membrane Properties. _J. Comput. Aided Mol. Des._ 2013, _27_ , 845–858
  * Pandit et al. 2003 Pandit, S. A.; Bostick, D.; Berkowitz, M. L. An Algorithm to Describe Molecular Scale Rugged Surfaces and its Application to the Study of a Water/Lipid Bilayer interface. _J. Chem. Phys._ 2003, _119_ , 2199–2205
  * King 1967 King, S. V. Ring Configurations in a Random Network Model of Vitreous Silica. _Nature_ 1967, _213_ , 1112–1113
  * Rahman and Stillinger 1973 Rahman, A.; Stillinger, F. H. Hydrogen-Bond Patterns in Liquid Water. _J. Am. Chem. Soc._ 1973, _95_ , 7943–7948
  * Guttman 1990 Guttman, L. Ring Structure of the Crystalline and Amorphous Forms of Silicon Dioxide. _J. Non-Cryst. Solids_ 1990, _116_ , 145–147
  * Franzblau 1991 Franzblau, D. S. Computation of Ring Statistics for Network Models of Solids. _Phys. Rev. B_ 1991, _44_ , 4925
  * Wooten 2002 Wooten, F. Structure, Odd Lines and Topological Entropy of Disorder of Amorphous Silicon. _Acta Cryst. A_ 2002, _58_ , 346–351
  * Yuan and Cormack 2002 Yuan, X.; Cormack, A. N. Efficient Algorithm for Primitive Ring Statistics in Topological Networks. _Comp. Mater. Sci._ 2002, _24_ , 343–360
  * Roux and Jund 2010 Roux, S. L.; Jund, P. Ring Statistics Analysis of topological Networks: New Approach and Application to Amorphous GeS2 and SiO2 Systems. _Comp. Mater. Sci._ 2010, _49_ , 70–83
  * Formanek and Martelli 2020 Formanek, M.; Martelli, F. Probing the Network Topology in Network–Forming Materials: the Case of Water. _AIP Adv._ 2020, _10_ , 055205
  * Luzar and Chandler 1996 Luzar, A.; Chandler, D. Hydrogen–Bond Kinetics in Liquid Water. _Nature_ 1996, _379_ , 55–57
  * Prada-Gracia et al. 2013 Prada-Gracia, D.; Shevchuk, R.; Rao, F. The Quest for Self–Consistency in Hydrogen Bond Definitions. _J. Chem. Phys._ 2013, _139_ , 084501
  * Shi et al. 2018 Shi, R.; Russo, J.; Tanaka, H. Common Microscopic Structural Origin for Water’s Thermodynamic and Dynamic Anomalies. _J. Chem. Phys._ 2018, _149_ , 224502
  * Goetzke and Klein 1991 Goetzke, K.; Klein, H.-J. Properties and Efficient Algorithmic Determination of Different Classes of Rings in Finite and Infinite Polyhedral Networks. _J. Non-Cryst. Solids_ 1991, _127_ , 215–220


6G

    sixth generation
AR

    auto-regressive
AoA

    angle of arrival
AoD

    angle of departure
AWGN

    additive white Gaussian noise
B5G

    beyond fifth generation
BS

    base station
BPSK

    binary phase-shift keying
IRS

    intelligent reflecting surface
mmWave

    millimeter Wave
MIMO

    multiple input multiple output
CSI

    channel state information
DFT

    discrete Fourier transform
HOSVD

    high order single value decomposition
LOS

    line of sight
MMSE

    minimum mean squared error
NMSE

    normalized mean square error
NLOS

    non-line of sight
SNR

    signal-to-noise ratio
SVD

    singular value decomposition
SE

    spectral efficiency
UE

    user equipment
ULA

    uniform linear array
URA

    uniform rectangular array
ALS

    alternating least squares
BALS

    bilinear alternating least squares
KF

    Kronecker factorization
KRF

    Khatri-Rao factorization
LS

    least squares
PARKRON

    PARAFAC-Khatri-Rao-Kronecker factorization
TBT

    Tucker-based tracking

# Tensor-based modeling/estimation of static channels in IRS-assisted MIMO
systems

Kenneth B. A. Benício, André L. F. de Almeida, Bruno Sokal, Fazal-E-Asim,
Behrooz Makki, and Gabor Fodor This work was supported by Ericsson Research,
Sweden, and Ericsson Innovation Center, Brazil, under UFC.51 Technical
Cooperation Contract Ericsson/UFC. This study was financed in part by
CAPES/Brazil - Finance Code 001, and CAPES/PRINT Proc. 88887.311965/2018-00.
André L. F. de Almeida thanks CNPq for its financial support under grant
312491/2020-4. G. Fodor was partially supported by the Digital Futures project
PERCy..

###### Abstract

This paper proposes a tensor-based parametric modeling and estimation
framework in MIMO systems assisted by intelligent reflecting surfaces (IRSs).
We present two algorithms that exploit the tensor structure of the received
pilot signal to estimate the concatenated channel. The first one is an
iterative solution based on the alternating least squares algorithm. In
contrast, the second method provides closed-form estimates of the involved
parameters using the high order single value decomposition. Our numerical
results show that our proposed tensor-based methods provide improved
performance compared to competing state-of-the-art channel estimation schemes,
thanks to the exploitation of the algebraic tensor structure of the combined
channel without additional computational complexity.

###### Index Terms:

channel estimation, intelligent reflecting surfaces, tensor-based algorithm,
complexity analysis

##  I Introduction

Over the last few years, intelligent reflecting surface (IRS) has been
considered one of the possible technologies to be deployed beyond fifth
generation (B5G) and sixth generation (6G) wireless networks due to its
potential to improve the coverage [1, 2]. An IRS is a $2$D panel composed of
many passive reflecting elements whose elements are capable of independently
changing the phase shifts of impinging electromagnetic waves to maximize the
signal-to-noise ratio (SNR) at the intended receiver [3]. Hence, channel
estimation must be performed at the end nodes of the network and the receiver
should estimate the involved channels from the received pilots reflected by
the IRS according to a training protocol. Several works have addressed this
problem, as mentioned in [4, 5, 6, 7, 8, 9, 10, 11].

As pointed out in [4], in IRS-assisted networks channel estimation methods can
be divided into structured and unstructured techniques exploiting the
parametric (geometric) modeling of the cascaded channel and methods exploiting
the combined channel structure, respectively. Also, in [5] a twin-IRS
structure consisting of two IRSs is proposed as a way to obtain the spatial
signatures of the involved channels. Also, authors in [6] propose a low-
complexity channel parameter estimation that exploits the decoupling of the
pilot design along the horizontal and vertical domains.

The authors in [7] propose channel parameter estimation using a low-rank
PARAFAC tensor in the context of millimeter Wave (mmWave) systems. The authors
in [8] use a tensor approach to perform supervised channel estimation, in
which the decoupling of the BS-IRS and IRS-UE channels is achieved. Then, [9]
proposes a tensor-based receiver formulated as a semi-blind problem that
jointly estimates the involved channels and transmitted data. The work in [10]
proposes a set of two tensor-based algorithms to do channel parameter
estimation under unknown IRS hardware impairments. In our previous work [11],
we propose a two-stage tensor-based framework for parametric channel parameter
estimation and data detection of time-varying channels based on a $4$th order
PARAFAC model.

In this paper, we propose a new signal modeling that exploits the geometric
channel structure to estimate the spatial signatures of the IRS-assisted
multiple input multiple output (MIMO) communication systems and formulate a
$3$rd order Tucker tensor model, and derive a set of two tensor algorithms
that solve the channel parameter estimation problem by either alternating
least squares (ALS) or high order single value decomposition (HOSVD).
Furthermore, we also study the computational complexity of the proposed
schemes and selected benchmark solutions. Our simulation results show that the
proposed techniques outperform the classic least squares (LS) and the state-
of-the-art Khatri-Rao factorization (KRF) [8] algorithms without increasing
the computational complexity.

Notation: Scalars, vectors, matrices, and tensors are represented as
$a,\bm{a},\bm{A}$, and $\mathcal{A}$. Also, $\bm{A}^{*}$, $\bm{A}^{\text{T}}$,
$\bm{A}^{\text{H}}$, and $\bm{A}^{\dagger}$ stand for the conjugate,
transpose, Hermitian, and pseudo-inverse, of a matrix $\bm{A}$, respectively.
The $j$th column of $\bm{A}\in\mathbb{C}^{I\times J}$ is denoted by
$\bm{a}_{j}\in\mathbb{C}^{I\times 1}$. The operator vec$(\cdot)$ transforms a
matrix into a vector by stacking its columns, e.g.,
$\text{vec}(\bm{A})=\bm{a}\in\mathbb{C}^{IJ\times 1}$, while the
unvec$(\cdot)_{{I\times J}}$ operator undo the operation. The operator
D$(\cdot)$ converts a vector into a diagonal matrix, $\text{D}_{j}(\bm{B})$
forms a diagonal matrix $R\times R$ out of the $j$th row of
$\bm{B}\in\mathbb{C}^{J\times R}$. Also, $\bm{I}_{N}$ denotes an identity
matrix of size $N\times N$. The symbols $\otimes$ and $\diamond$ indicate the
Kronecker and Khatri-Rao products.

##  II System Model

Figure 1: Proposed IRS-assisted MIMO system scenario.

We consider an uplink IRS-assisted MIMO scenario with a base station (BS)
equipped with $M$ receiver antennas, which receives a signal from a user
equipment (UE) equipped with $Q$ transmit antennas via a passive IRS with $N$
reflecting elements as shown in Fig. 1. The transmission has length of $T$
time-slots and the received pilot signal at the $t$-th time-slot is given by

$\displaystyle\bm{y}_{t}$
$\displaystyle=\bm{G}\text{D}(\bm{s}_{t})\bm{H}\bm{z}_{t}+\bm{v}_{t}\in\mathbb{C}^{M\times
1},$ (1)

where $\bm{z}_{t}$ is the pilot sequence, $\text{D}(\bm{s}_{t})$ is the IRS
phase-shift matrix, and $\bm{v}_{t}$ is the additive white Gaussian noise
(AWGN) vector with $t\in\\{1,\cdots,T\\}$. We assume that the IRS-UE channel,
$\bm{H}$, and the BS-IRS channel, $\bm{G}$, remain constant during $T$ time-
slots and consider a mmWave scenario adopting a multipath channel model [12]
for the involved channels. We can express these channel matrices as follows

$\displaystyle\bm{G}$
$\displaystyle=\sqrt{\frac{K_{G}}{K_{G}+1}}\bm{G}^{(\text{LOS})}+\sqrt{\frac{1}{K_{G}+1}}\bm{G}^{\text{(NLOS)}},$
$\displaystyle\bm{H}$
$\displaystyle=\sqrt{\frac{K_{H}}{K_{H}+1}}\bm{H}^{(\text{LOS})}+\sqrt{\frac{1}{K_{H}+1}}\bm{H}^{\text{(NLOS)}},$

with

$\displaystyle\bm{G}^{(\text{LOS})}$
$\displaystyle=\alpha^{(1)}\bm{a}_{\text{r}_{\text{x}}}(\mu^{(1)}_{\text{bs}})\bm{b}^{\text{(irs)}\text{H}}_{\text{t}_{\text{x}}}(\mu^{(1)}_{\text{irs}_{{\text{D}}}},\psi^{(1)}_{\text{irs}_{\text{D}}}),$
(2) $\displaystyle\bm{G}^{(\text{NLOS})}$
$\displaystyle=\sum^{L_{1}}_{l_{1}=2}\alpha^{(l_{1})}\bm{a}_{\text{r}_{\text{x}}}(\mu^{(l_{1})}_{\text{bs}})\bm{b}^{\text{(irs)}\text{H}}_{\text{t}_{\text{x}}}(\mu^{(l_{1})}_{\text{irs}_{{\text{D}}}},\psi^{(l_{1})}_{\text{irs}_{\text{D}}}),$
(3) $\displaystyle\bm{H}^{(\text{LOS})}$
$\displaystyle=\beta^{(1)}\bm{b}^{\text{(irs)}}_{\text{r}_{\text{x}}}(\mu^{(1)}_{\text{irs}_{{\text{A}}}},\psi^{(1)}_{\text{irs}_{\text{A}}})\bm{a}_{\text{t}_{\text{x}}}^{\text{H}}(\mu^{(1)}_{\text{ue}}),$
(4) $\displaystyle\bm{H}^{(\text{NLOS})}$
$\displaystyle=\sum^{L_{2}}_{l_{2}=2}\beta^{(l_{2})}\bm{b}^{\text{(irs)}}_{\text{r}_{\text{x}}}(\mu^{(l_{2})}_{\text{irs}_{{\text{A}}}},\psi^{(l_{2})}_{\text{irs}_{\text{A}}})\bm{a}_{\text{t}_{\text{x}}}^{\text{H}}(\mu^{(l_{2})}_{\text{ue}}),$
(5)

where $\bm{G}^{(\text{LOS})}$, and $\bm{H}^{(\text{LOS})}$, are the line of
sight (LOS) components and $\bm{G}^{(\text{NLOS})}$, and
$\bm{H}^{(\text{NLOS})}$ are the non-line of sight (NLOS) components. Also,
$K_{G}$ and $K_{H}$ are the Rician factors for channels $\bm{G}$ and $\bm{H}$.
The $l$th one-dimensional steering vector of the BS is
$\bm{a}_{\text{r}_{\text{x}}}(\mu^{(l_{1})}_{\text{bs}})$ having spatial
frequency defined as
$\mu^{(l_{1})}_{\text{bs}}=\pi\text{cos}(\phi^{(l_{1})}_{\text{bs}})$ with
$\phi^{(l_{1})}_{\text{bs}}$ being the angle of arrival (AoA), which can be
further written in terms of spatial frequency assuming uniform linear array
(ULA) as [13]

$\displaystyle\bm{a}_{\text{r}_{\text{x}}}(\mu^{(l_{1})}_{\text{bs}})=\left[1,\cdots,e^{-j\pi(M-1)\mu^{(l_{1})}_{\text{bs}}}\right]^{\text{T}}\in\mathbb{C}^{M\times
1}.$ (6)

Similarly, the $p$th one-dimensional steering vector for the UE is
$\bm{a}_{\text{t}_{\text{x}}}(\mu^{(l_{2})}_{\text{ue}})$ having spatial
frequency, which is defined as
$\mu^{(l_{2})}_{\text{ue}}=\pi\text{cos}(\phi^{(l_{2})}_{\text{ue}})$, with
$\phi^{(l_{2})}_{\text{ue}}$ being the angle of departure (AoD), and can be
written in terms of spatial frequency as

$\displaystyle\bm{a}_{\text{t}_{\text{x}}}(\mu^{(l_{2})}_{\text{ue}})=\left[1,\cdots,e^{-j\pi(Q-1)\mu^{(l_{2})}_{\text{ue}}}\right]^{\text{T}}\in\mathbb{C}^{Q\times
1}.$ (7)

At the IRS,
$\bm{b}^{(\text{irs})}_{\text{r}_{\text{x}}}(\mu^{(l_{2})}_{\text{irs}_{{\text{A}}}},\psi^{(l_{2})}_{\text{irs}_{\text{A}}})$
is the 2D steering vector having spatial frequencies defined as
$\mu^{(l_{2})}_{\text{irs}_{{\text{A}}}}=\pi\text{cos}(\phi^{(l_{2})}_{\text{irs}_{\text{A}}})\text{sin}(\theta^{(l_{2})}_{\text{irs}_{\text{A}}})$
and
$\psi^{(l_{2})}_{\text{irs}_{\text{A}}}=\pi\text{cos}(\phi^{(l_{2})}_{\text{irs}_{\text{A}}})$,
where $\phi^{(l_{2})}_{\text{irs}_{\text{A}}}$ and
$\theta^{(l_{2})}_{\text{irs}_{\text{A}}}$ are the azimuth AoA and the
elevation AoA, respectively. This can be further written as the Kronecker
product between two steering vectors as [13]

$\displaystyle\bm{b}^{\text{(irs)}}_{\text{r}_{\text{x}}}(\mu^{(l_{2})}_{\text{irs}_{{\text{A}}}},\psi^{(l_{2})}_{\text{irs}_{\text{A}}})=\bm{b}^{\text{(irs)}}_{\text{r}_{\text{x}}}(\mu^{(l_{2})}_{\text{irs}_{{\text{A}}}})\otimes\bm{b}^{\text{(irs)}}_{\text{r}_{\text{x}}}(\psi^{(l_{2})}_{\text{irs}_{\text{A}}})\in\mathbb{C}^{N\times
1}.$ (8)

The IRS transmission steering vector,
$\bm{b}^{\text{(irs)\text{H}}}_{\text{t}_{\text{x}}}(\mu^{(l_{1})}_{\text{irs}_{{\text{D}}}},\psi^{(l_{1})}_{\text{irs}_{\text{D}}})$,
is defined similarly. The IRS phase-shift vector is defined as
$\bm{s}_{t}=\left[e^{j\theta_{t}},\cdots,e^{j\theta_{N,t}}\right]^{\text{T}}\in\mathbb{C}^{N\times
1}$, where $\theta_{n,t}$ is the phase-shift of the $n$th IRS element at the
$t$th time slot. Also,
$\bm{\alpha}=\left[\alpha^{(1)},\cdots,\alpha^{(L_{1})}\right]^{\text{T}}\in\mathbb{C}^{L_{1}\times
1}$ and
$\bm{\beta}=\left[\beta^{(1)},\cdots,\beta^{(L_{2})}\right]^{\text{T}}\in\mathbb{C}^{L_{2}\times
1}$ represent the path loss and fading components of the BS-IRS and IRS-UE
channels, respectively. Both channels are written in matrix notation as

$\displaystyle\bm{G}$
$\displaystyle=\bm{A}_{\text{r}_{\text{x}}}\bm{D}(\bm{\alpha})\bm{B}^{\text{H}}_{\text{t}_{\text{x}}}\in\mathbb{C}^{M\times
N},$ (9) $\displaystyle\bm{H}$
$\displaystyle=\bm{B}_{\text{r}_{\text{x}}}\bm{D}(\bm{\beta})\bm{A}_{\text{t}_{\text{x}}}^{\text{H}}\in\mathbb{C}^{N\times
Q},$ (10)

where $\bm{A}_{\text{r}_{\text{x}}}$ and $\bm{B}_{\text{r}_{\text{x}}}$ are
the steering matrices defined as

$\displaystyle\bm{A}_{\text{r}_{\text{x}}}$
$\displaystyle=\left[\bm{a}_{\text{r}_{\text{x}}}(\mu^{(1)}_{\text{bs}}),\cdots,\bm{a}_{\text{r}_{\text{x}}}(\mu^{(L_{1})}_{\text{bs}})\right]\in\mathbb{C}^{M\times
L_{1}},$ $\displaystyle\bm{B}_{\text{r}_{\text{x}}}$
$\displaystyle=\left[\bm{b}^{(\text{irs})}_{\text{r}_{\text{x}}}(\mu^{(1)}_{\text{irs}_{{\text{A}}}},\psi^{(1)}_{\text{irs}_{\text{A}}}),\cdots,\bm{b}^{\text{(irs)}}_{\text{r}_{\text{x}}}(\mu^{(L_{2})}_{\text{irs}_{{\text{A}}}},\psi^{(L_{2})}_{\text{irs}_{\text{A}}})\right]\in\mathbb{C}^{N\times
L_{2}},$

with $\bm{B}_{\text{t}_{\text{x}}}$ and $\bm{A}_{\text{t}_{\text{x}}}$ being
defined in similar manner.

##  III Pilot-based Parameter Estimation

In this section, we describe the proposed tensor-based methods for channel
parameter estimation, namely Tucker-ALS as in Alg. 1, and Tucker-HOSVD as in
Alg.2. The main idea is to exploit the geometric structure of the involved
channels by using a tensor approach.

###  III-A Tensor-Based Parameter Estimation

In this section, we formulate a tensor-based approach to estimate the channel
parameters. Using
$\text{vec}(\bm{A}\bm{B}\bm{C})=(\bm{C}^{\text{T}}\otimes\bm{A})\text{vec}(\bm{B})$
and
$\text{vec}(\bm{A}\text{D}\left(\bm{b}\right)\bm{C})=(\bm{C}^{\text{T}}\diamond\bm{A})\bm{b}$
in (1), yields

$\displaystyle\bm{y}_{t}$ $\displaystyle=\text{vec
}(\bm{I}_{M}\bm{G}\text{D}(\bm{s}_{t})\bm{H}\bm{z}_{t})+\bm{v}_{t}\in\mathbb{C}^{M\times
1},$
$\displaystyle=(\bm{s}^{\text{T}}_{t}\otimes\bm{z}^{\text{T}}_{t}\otimes\bm{I}_{M})\text{vec}(\bm{H}^{\text{T}}\diamond\bm{G})+\bm{v}_{t}.$
$\displaystyle\bm{y}_{t}$ $\displaystyle=\text{vec
}(\bm{I}_{M}\bm{G}\text{D}(\bm{s}_{t})\bm{H}\bm{z}_{t})+\bm{v}_{t}\in\mathbb{C}^{M\times
1},$
$\displaystyle=(\bm{z}^{\text{T}}_{t}\otimes\bm{I}_{M})\text{vec}(\bm{G}\text{D}(\bm{s}_{t})\bm{H})+\bm{v}_{t},$
$\displaystyle=(\bm{z^{\text{T}}}_{t}\otimes\bm{I}_{M})(\bm{H}^{\text{T}}\diamond\bm{G})\bm{s}_{t}+\bm{v}_{t},$

and, applying the first property again, we have

$\displaystyle\bm{y}_{t}$
$\displaystyle=\text{vec}[(\bm{z}^{\text{T}}_{t}\otimes\bm{I}_{M})(\bm{H}^{\text{T}}\diamond\bm{G})\bm{s}_{t}]+\bm{v}_{t},$
$\displaystyle=(\bm{s}^{\text{T}}_{t}\otimes\bm{z}^{\text{T}}_{t}\otimes\bm{I}_{M})\text{vec}(\bm{H}^{\text{T}}\diamond\bm{G})+\bm{v}_{t}.$

Collecting the signals during the $T$ symbol periods yields

$\displaystyle\bm{y}$
$\displaystyle=\left[\bm{y}_{1}^{\text{T}},\cdots,\bm{y}_{T}^{\text{T}}\right]^{\text{T}},$
$\displaystyle=[(\bm{S}\diamond\bm{Z})^{\text{T}}\otimes\bm{I}_{M}]\text{vec}(\bm{H}^{\text{T}}\diamond\bm{G})+\bm{v},$
$\displaystyle=\bm{\Omega}\hskip 4.26773pt\bm{u}+\bm{v}\in\mathbb{C}^{MT\times
1},$ (11)

where $\bm{S}=\left[\bm{s}_{1},\cdots,\bm{s}_{T}\right]\in\mathbb{C}^{N\times
T}$, $\bm{Z}=\left[\bm{z}_{1},\cdots,\bm{z}_{T}\right]\in\mathbb{C}^{Q\times
T}$ are matrices collecting the IRS phase-shifts and pilots,
$\bm{\Omega}=(\bm{S}\diamond\bm{Z})^{\text{T}}\otimes\bm{I}_{M}\in\mathbb{C}^{MT\times
MQN}$,
$\bm{u}=\textrm{vec}(\bm{H}^{\text{T}}\diamond\bm{G})\in\mathbb{C}^{MQN\times
1}$, and
$\bm{v}=\left[\bm{v}^{\text{T}}_{1},\cdots,\bm{v}^{\text{T}}_{T}\right]^{\text{T}}\in\mathbb{C}^{MT\times
1}$ is the AWGN noise term. From (11), we obtain the following LS problem

$\displaystyle\bm{\hat{u}}=\underset{\bm{u}}{\text{arg
min}}\left|\left|\bm{y}-\bm{\Omega}\bm{u}\right|\right|^{2}_{2},$ (12)

where the solution requires $T\geq QN$ and is given by

$\displaystyle\bm{\hat{u}}$
$\displaystyle=\bm{\Omega}^{\dagger}\bm{y}\in\mathbb{C}^{MQN\times 1}.$ (13)

Let us define $\bm{R}=\textrm{unvec}_{MQ\times
N}(\bm{\hat{u}})\approx\bm{H}^{\text{T}}\diamond\bm{G}\in\mathbb{C}^{MQ\times
N}$, where the approximation is exact in a noiseless scenario. Using (9) and
(10), while applying property
$(\bm{AC})\diamond(\bm{BD})=(\bm{A}\otimes\bm{B})(\bm{C}\diamond\bm{D})$, we
have

$\displaystyle\bm{R}$
$\displaystyle\approx[\bm{A}_{\text{t}_{\text{x}}}^{*}\text{D}(\bm{\beta})\bm{B}^{\text{T}}_{\text{r}_{\text{x}}}]\diamond[\bm{A}_{\text{r}_{\text{x}}}\text{D}(\bm{\alpha})\bm{B}^{\text{H}}_{\text{t}_{\text{x}}}],$
$\displaystyle\approx(\bm{A}_{\text{t}_{\text{x}}}^{*}\otimes\bm{A}_{\text{r}_{\text{x}}})[\text{D}(\bm{\beta})\otimes\text{D}(\bm{\alpha})](\bm{B}^{\text{T}}_{\text{r}_{\text{x}}}\diamond\bm{B}^{\text{H}}_{\text{t}_{\text{x}}}).$
(14)

Defining $\bm{f}=\bm{\beta}\otimes\bm{\alpha}\in\mathbb{C}^{L_{1}L_{2}\times
1}$, (14) can be expressed as

$\displaystyle\bm{R}$
$\displaystyle\approx(\bm{A}_{\text{t}_{\text{x}}}^{*}\otimes\bm{A}_{\text{r}_{\text{x}}})\text{D}(\bm{f})\bm{P}^{\text{T}}_{B}\in\mathbb{C}^{MQ\times
N},$ (15)

where
$\bm{P}_{B}=(\bm{B}^{\text{T}}_{\text{r}_{\text{x}}}\diamond\bm{B}^{\text{H}}_{\text{t}_{\text{x}}})\in\mathbb{C}^{N\times
L_{1}L_{2}}$ is the IRS geometry information. Note that
$\text{D}(\bm{f})\in\mathbb{C}^{L_{1}L_{2}\times L_{1}L_{2}}$ can be viewed as
the $3$-mode unfolding of the tensor $\mathcal{F}\in\mathbb{C}^{L_{1}\times
L_{2}\times L_{1}L_{2}}$, i.e.,
$\text{D}(\bm{f})=\left[\mathcal{F}\right]_{(3)}$. The tensor $\mathcal{F}$ is
given by [14]

$\displaystyle\mathcal{F}$
$\displaystyle=\left(\mathcal{I}_{3,L_{2}}\otimes^{2,3}_{2,3}\mathcal{I}_{3,L_{1}}\right)\times_{3}\bm{f}^{\text{T}}\in\mathbb{C}^{L_{1}\times
L_{2}\times L_{1}L_{2}},$ (16)

where $\mathcal{I}_{3,L_{2}}$ and $\mathcal{I}_{3,L_{1}}$ are identity tensors
and $\otimes^{2,3}_{2,3}$ is the selective Kronecker product (SKP) [14], from
which we arranged (15) following a third-way Tucker tensor structure
$\mathcal{R}\in\mathbb{C}^{M\times Q\times N}$, as shown in Fig. 2, in terms
of $n$-mode product as

$\displaystyle\mathcal{R}$
$\displaystyle\approx\mathcal{F}\times_{1}\bm{A}_{\text{r}_{\text{x}}}\times_{2}\bm{A}^{*}_{\text{t}_{\text{x}}}\times_{3}\bm{P}_{B},$
(17)

with the matrix unfoldings of $\mathcal{R}$ given by

$\displaystyle\left[\mathcal{R}\right]_{(1)}$
$\displaystyle\approx\bm{A}_{\text{r}_{\text{x}}}\left[\mathcal{F}\right]_{(1)}(\bm{P}_{B}\otimes\bm{A}^{*}_{\text{t}_{\text{x}}})^{\text{T}}\in\mathbb{C}^{M\times
QN},$ $\displaystyle\left[\mathcal{R}\right]_{(2)}$
$\displaystyle\approx\bm{A}^{*}_{\text{t}_{\text{x}}}\left[\mathcal{F}\right]_{(2)}(\bm{P}_{B}\otimes\bm{A}_{\text{r}_{\text{x}}})^{\text{T}}\in\mathbb{C}^{Q\times
MN},$ $\displaystyle\left[\mathcal{R}\right]_{(3)}$
$\displaystyle\approx\bm{P}_{B}\left[\mathcal{F}\right]_{(3)}(\bm{A}^{*}_{\text{t}_{\text{x}}}\otimes\bm{A}_{\text{r}_{\text{x}}})^{\text{T}}\in\mathbb{C}^{N\times
MQ}.$

Consequently, the estimation of $\bm{A}_{\text{r}_{\text{x}}}$,
$\bm{A}_{\text{t}_{\text{x}}}$, $\bm{P}_{B}$ and $\mathcal{F}$ consists of
solving the following problem

$\left\\{\hat{\bm{A}}_{\text{r}_{\text{x}}},\hat{\bm{A}}_{\text{t}_{\text{x}}},\hat{\bm{P}_{B}},\hat{\mathcal{F}}\right\\}=\underset{\bm{A}_{\text{r}_{\text{x}}},\bm{A}_{\text{t}_{\text{x}}},\bm{P}_{B},\mathcal{F}}{\text{arg
min}}\left|\left|\begin{split}\mathcal{R}-\mathcal{F}\times_{1}\bm{A}_{\text{r}_{\text{x}}}\\\
\times_{2}\bm{A}^{*}_{\text{t}_{\text{x}}}\times_{3}\bm{P}_{B}\end{split}\right|\right|^{2}_{\text{F}},$
(18)

which can be performed by means of the well-known ALS algorithm [15] or the
HOSVD algorithm [16].

Figure 2: Proposed Tucker-based model.

###  III-B Tucker ALS

From (17) we can derive an iterative solution based on the well-known ALS
algorithm [15]. Here, the algorithm consists of an estimate
$\bm{A}_{\text{r}_{\text{x}}}$, $\bm{A}_{\text{t}_{\text{x}}}$, $\bm{P}_{B}$,
and $\mathcal{F}$ in an alternating way by iteratively solving the following
cost functions

$\displaystyle\hat{\bm{A}}_{\text{r}_{\text{x}}}$
$\displaystyle=\underset{\bm{A}_{\text{r}_{\text{x}}}}{\text{arg
min}}\left|\left|\left[\mathcal{R}\right]_{(1)}-\bm{A}_{\text{r}_{\text{x}}}\left[\mathcal{F}\right]_{(1)}(\bm{P}_{B}\otimes\bm{A}^{*}_{\text{t}_{\text{x}}})^{\text{T}}\right|\right|^{2}_{\text{F}},$
(19) $\displaystyle\hat{\bm{A}}_{\text{t}_{\text{x}}}$
$\displaystyle=\underset{\bm{A}_{\text{t}_{\text{x}}}}{\text{arg
min}}\left|\left|\left[\mathcal{R}\right]_{(2)}-\bm{A}^{*}_{\text{t}_{\text{x}}}\left[\mathcal{F}\right]_{(2)}(\bm{P}_{B}\otimes\bm{A}_{\text{r}_{\text{x}}})^{\text{T}}\right|\right|^{2}_{\text{F}},$
(20) $\displaystyle\hat{\bm{P}}_{B}$
$\displaystyle=\underset{\bm{P}_{B}}{\text{arg
min}}\left|\left|\left[\mathcal{R}\right]_{(3)}-\bm{P}_{B}\left[\mathcal{F}\right]_{(3)}(\bm{A}^{*}_{\text{t}_{\text{x}}}\otimes\bm{A}_{\text{r}_{\text{x}}})^{\text{T}}\right|\right|^{2}_{\text{F}},$
(21) $\displaystyle\hat{\bm{f}}$ $\displaystyle=\underset{\bm{f}}{\text{arg
min}}\left|\left|\text{vec}(\left[\mathcal{R}\right]_{(3)})-[(\bm{A}^{*}_{\text{t}_{\text{x}}}\otimes\bm{A}_{\text{r}_{\text{x}}})\diamond\bm{P}_{B}]\bm{f}\right|\right|^{2}_{2},$
(22)

where $\bm{f}$ is related to the diagonal elements of
$\left[\mathcal{F}\right]_{(3)}$ and its expression is derived by applying
$\text{vec}(\bm{A}\text{D}\left(\bm{b}\right)\bm{C})=(\bm{C}^{\text{T}}\diamond\bm{A})\bm{b}$
to $\left[\mathcal{R}\right]_{(3)}$. The solutions for (19)-(22) are
respectively given by

$\displaystyle\hat{\bm{A}}_{\text{r}_{\text{x}}}$
$\displaystyle=\left[\mathcal{R}\right]_{(1)}\left[\left[\mathcal{F}\right]_{(1)}(\bm{P}_{B}\otimes\bm{A}^{*}_{\text{t}_{\text{x}}})^{\text{T}}\right]^{\dagger},$
(23) $\displaystyle\hat{\bm{A}}_{\text{t}_{\text{x}}}$
$\displaystyle=\left[\mathcal{R}\right]_{(2)}\left[\left[\mathcal{F}\right]_{(2)}(\bm{P}_{B}\otimes\bm{A}_{\text{r}_{\text{x}}})^{\text{T}}\right]^{\dagger},$
(24) $\displaystyle\hat{\bm{P}}_{B}$
$\displaystyle=\left[\mathcal{R}\right]_{(3)}\left[\left[\mathcal{F}\right]_{(3)}(\bm{A}^{*}_{\text{t}_{\text{x}}}\otimes\bm{A}_{\text{r}_{\text{x}}})^{\text{T}}\right]^{\dagger},$
(25) $\displaystyle\hat{\bm{f}}$
$\displaystyle=\left[(\bm{A}^{*}_{\text{t}_{\text{x}}}\otimes\bm{A}_{\text{r}_{\text{x}}})\diamond\bm{P}_{B}\right]^{\dagger}\text{vec}(\left[\mathcal{R}\right]_{(3)}),$
(26)

with each solution requiring that $QN\geq L_{1}$, $MN\geq L_{2}$, $MQ\geq
L_{1}L_{2}$, and $NMQ\geq L_{1}L_{2}$. These four conditions are related to
the convergence of the LS estimates $\bm{A}_{\text{r}_{\text{x}}}$,
$\bm{A}^{*}_{\text{t}_{\text{x}}}$, $\bm{P}_{B}$, and $\bm{f}$, respectively.
The proposed ALS algorithm consists of four iterative and alternating update
steps that follow the LS solutions (23)-(26). At each update, the
reconstruction error is minimized according to one given factor matrix by
fixing the other matrices to their estimation obtained at the previous update.
This procedure is repeated until the convergence is acknowledged, which
happens when the reconstruction error, given by
$e(i)=||\mathcal{R}-\hat{\mathcal{R}}(i)||^{2}_{\text{F}}$, achieves
$||e(i)-e(i-1)||\leq\epsilon$ and $\epsilon$ is the threshold parameter with
$\hat{\mathcal{R}}(i)$ being the reconstruct tensor model at the $i$th
iteration. In this work, we initialize the factor matrices randomly and the
convergence threshold is set to $\epsilon=10^{-5}$. It is worth noting that
our solution using the ALS is not unique although, after the convergence, the
intrinsic scaling and permutation ambiguities disappear on the estimated
tensor model
$\hat{\mathcal{R}}=\hat{\mathcal{F}}\times_{1}\hat{\bm{A}}_{\text{r}_{\text{x}}}\times_{2}\hat{\bm{A}}^{*}_{\text{t}_{\text{x}}}\times_{3}\hat{\bm{P}}_{B}$.

Algorithm 1 ALS

1:Tensor $\mathcal{R}$

2:while $||e(i)-e(i-1)||\geq\delta$ do

3: Find a LS estimate of $\bm{A}_{\text{r}_{\text{x}}}$ with (23).

4: Find a LS estimate of $\bm{A}_{\text{t}_{\text{x}}}$ with (24).

5: Find a LS estimate of $\bm{P}_{B}$ with (25).

6: Find an estimate of $\bm{f}$ with (26).

7: Build $\mathcal{F}$ with (16).

8: Repeat until convergence.

9:end while

10:return
$\hat{\mathcal{R}}=\hat{\mathcal{F}}\times_{1}\hat{\bm{A}}_{\text{r}_{\text{x}}}\times_{2}\hat{\bm{A}}^{*}_{\text{t}_{\text{x}}}\times_{3}\hat{\bm{P}}_{B}$

###  III-C Tucker HOSVD

Considering the Tucker model in (17), we can also estimate its factors by
finding a multi-linear rank ($L_{1},L_{2},L_{1}L_{2}$) approximation to
$\mathcal{R}$. This can be done by means of the state-of-art truncated HOSVD
algorithm. This algorithm consists of computing multiple singular value
decompositions (SVDs), one for each unfolding of $\mathcal{R}$, as follows

$\displaystyle\left[\mathcal{R}\right]_{(1)}$
$\displaystyle=\bm{U}^{(1)}\bm{\Sigma}^{(1)}\bm{V}^{(1)\text{H}},$ (27)
$\displaystyle\left[\mathcal{R}\right]_{(2)}$
$\displaystyle=\bm{U}^{(2)}\bm{\Sigma}^{(2)}\bm{V}^{(2)\text{H}},$ (28)
$\displaystyle\left[\mathcal{R}\right]_{(3)}$
$\displaystyle=\bm{U}^{(3)}\bm{\Sigma}^{(3)}\bm{V}^{(3)\text{H}}.$ (29)

The estimates of the steering matrices $\bm{A}_{\text{r}_{\text{x}}}$,
$\bm{A}_{\text{t}_{\text{x}}}$, and $\bm{P}_{B}$ are found from the dominant
$L_{1}$, $L_{2}$, and $L_{1}L_{2}$ left singular vectors of
$\left[\mathcal{R}\right]_{(1)}$, $\left[\mathcal{R}\right]_{(2)}$, and
$\left[\mathcal{R}\right]_{(3)}$, respectively

$\displaystyle\hat{\bm{A}}_{\text{r}_{\text{x}}}=\bm{U}^{(1)}_{.1:L_{1}},\hat{\bm{A}}_{\text{t}_{\text{x}}}=\bm{U}^{(2)}_{.1:L_{2}},\hat{\bm{P}}_{B}=\bm{U}^{(3)}_{.1:L_{1}L_{2}}.$
(30)

Finally, an estimate of the core tensor in (17) is obtained as

$\displaystyle\hat{\bm{f}}$
$\displaystyle=\left[(\hat{\bm{A}}^{*}_{\text{t}_{\text{x}}}\otimes\hat{\bm{A}}_{\text{r}_{\text{x}}})\diamond\hat{\bm{P}}_{B}\right]^{\dagger}\text{vec}(\left[\mathcal{R}\right]_{(3)}),$
(31)

where $\bm{f}$. Similar to the ALS solution, the HOSVD is not unique once any
transformation applied to the core tensor does not change the tensor fit [16]
except in the case where the core tensor is known [17].

Algorithm 2 Tucker HOSVD

1:Tensor $\mathcal{R}$

2:Define the singular value decomposition (SVD)s of
$\left[\mathcal{R}\right]_{(1)}$, $\left[\mathcal{R}\right]_{(2)}$, and
$\left[\mathcal{R}\right]_{(3)}$ with (27)-(29).

3:Compute the estimation of matrices $\bm{A}_{\text{r}_{\text{x}}}$,
$\bm{A}_{\text{t}_{\text{x}}}$, and $\bm{P}_{B}$ with (30).

4:Compute the estimation of $\bm{f}=\text{vec}(\mathcal{F})$ with (31)

5:Reconstruct $\hat{\mathcal{F}}$ with (16).

6:return
$\hat{\mathcal{R}}=\hat{\mathcal{F}}\times_{1}\hat{\bm{A}}_{\text{r}_{\text{x}}}\times_{2}\hat{\bm{A}}^{*}_{\text{t}_{\text{x}}}\times_{3}\hat{\bm{P}}_{B}$

###  III-D Computational Complexity

In Table LABEL:tab:complexity, we describe the computational complexity for
the selected benchmark algorithms, LS (13) and KRF [8], and our proposed
solutions, ALS and HOSVD. Consider that the pseudo-inverse of a matrix
$\bm{A}\in\mathbb{C}^{I\times J}$, with $I>J$, and its rank-$R$ SVD
approximation have complexities $\mathcal{O}(IJ^{2})$ and $\mathcal{O}(IJR)$,
respectively. Since the design of $\bm{\Omega}$ at (13) is orthogonal, we have
$\bm{\Omega}^{\dagger}=\bm{\Omega}^{\text{H}}$ which lower the cost of the
pseudo-inverse computation. The KRF [8] estimates the combined channel,
$\bm{R}=\bm{H}^{\text{T}}\diamond\bm{G}$, by finding estimates of both
$\bm{G}$ and $\bm{H}$ that solves a set of $N$ rank-one approximations using
the SVD. Regarding the proposed algorithms, the ALS computes $4$ pseudo-
inverses (23)-(26) along $\text{\acs{ALS}}_{\text{iter}}$ iterations until
convergence, while the HOSVD involves $3$ SVDs (27)-(29) and a pseudo-inverse
(31).

##  IV Simulation Results

We evaluate the performance of the proposed tensor-based algorithm by
comparing it with the reference parameter estimation method based on the KRF
[8]. The pilot matrix $\bm{Z}\in\mathbb{C}^{Q\times T}$ is designed as a
Hadamard matrix, while a discrete Fourier transform (DFT) is adopted for the
IRS phase-shift matrix $\bm{S}$. The angular parameters
$\phi^{(l_{1})}_{\text{bs}}$ and $\phi^{(l_{2})}_{\text{ue}}$ are randomly
generated from a uniform distribution between $[-\pi,\pi]$ while the IRS
elevation and azimuth angles of arrival and departure are randomly generated
from a uniform distribution between $[-\pi/2,\pi/2]$. The fading coefficients
$\bm{\alpha}$ and $\bm{\beta}$ are modeled as independent Gaussian random
variables $\mathcal{CN}(0,1)$. The parameter estimation accuracy is evaluated
in terms of the normalized mean square error (NMSE) given by

$\displaystyle\text{\acs{NMSE}}(\bm{R})=\mathbb{E}\left\\{\frac{\left|\left|\bm{R}^{(m)}-\hat{\bm{R}}^{(m)}\right|\right|^{2}_{\text{F}}}{\left|\left|\bm{R}^{(m)}\right|\right|^{2}_{\text{F}}}\right\\},$

with $\bm{R}=(\bm{H}^{\text{T}}\diamond\bm{G})$ being the estimated channel at
the $m$th experiment, $M=10^{4}$ being the number of Monte Carlo experiments,
and $\sigma^{2}_{\bm{V}}$ is the noise variance. Unless otherwise stated, the
training SNR is $30$ dB, the Rician factor of the LOS channel is $K_{G}=10$
dB, and the Rician factor of the NLOS channel is $K_{H}=-10$ dB. At Figs. 4,
4, and 7 we assume $\\{M=4,Q=4,L_{1}=1,L_{2}=4,N=16,\text{ and }T=64\\}$, at
Fig. 6 we assume $\\{M=8,Q=8,N=16,\text{ and }T=128\\}$, and finally at Fig. 6
we assume $\\{M=8,Q=8,L_{1}=2,L_{2}=2,\text{ and }T=8N\\}$.

In Fig. 4, we show the impact of the Rician factors, $K_{G}$ and $K_{H}$, on
the NMSE performance associated with the estimation of the combined channel
$\bm{R}=\bm{H}^{\text{T}}\diamond\bm{G}$, for the proposed algorithms and the
benchmark KRF [8]. The performance improvements with increasing values of the
Rician factor come from the fact that as $K_{G}$ and $K_{H}$ increases, both
channels are dominated by a LOS component, which happens because the channel
approaches a unitary rank and the noise rejection coming from the SVD. We can
also see that the proposed algorithms outperform the classic LS (13) and the
state-of-the-art KRF [8] algorithms by approximately $10$ dB and $5$ dB,
respectively. Moreover, the ALS performance is better than that of the HOSVD
in about 1 dB. In Fig. 4, we evaluate the NMSE performance as a function of
the training SNR and verify the same gains as the evaluation of the NMSE in
terms of the Rician factors. In all considered methods, we observe that the
NMSE decreases with SNR linearly.

In Fig. 6, we evaluate the number of iterations required by the proposed ALS
algorithm to converge as a function of the SNR for a varying number of paths,
$L_{1}$ and $L_{2}$. In this figure, we define $\epsilon=10^{-5}$ as the
target convergence criterion, meaning that the algorithm convergence is
declared when the error between consecutive iterations is less than
$\epsilon$. As expected, in the low SNR region the ALS algorithm takes more
iterations to achieve convergence as the total number of components, i.e., the
product $L_{1}L_{2}$, increases. At the high SNR region, the required number
of iterations is the same. In Fig. 6, we observe that, as the number of
reflecting elements $N$ increases, fewer iterations are needed for
convergence. This is linked to LS (13) since, if $N$ increases, we sense the
channel longer because of condition $T\geq QN$.

In Fig. 7, we analyze the computational complexity of the benchmark
algorithms, LS (13) and KRF [8], and the proposed algorithms, ALS and HOSVD,
for fixed parameters $M,Q,L_{1},L_{2}$ while varying $N$ according to Table
LABEL:tab:complexity. To compute the cost of the KRF [8], the ALS, and the
HOSVD, we take into account the extra cost of the LS (13) step, which is the
most complex operation of the competing and our proposed solutions (see Table
LABEL:tab:complexity). We observe that the competing LS and KRF [8] algorithms
have approximately the same cost as the proposed ALS and HOSVD solutions.

Figure 3: NMSE as a function of the Rician factors $K_{G}$ and $K_{H}$.

Figure 4: NMSE as a function of the training SNR.

Figure 5: Iterations to converge as function of the SNR and for increasing
$L_{1}$ and $L_{2}$.

Figure 6: Iterations to converge as function of the SNR and for increasing
$N$.

Figure 7: Computational complexity of the proposed solutions and the benchmark
algorithms LS (13) KRF [8].

##  V Conclusions

This paper proposes two tensor-based channel parameter estimation algorithms
in IRS-aided MIMO communications. From our simulation results, we observed
that proposed ALS and HOSVD algorithms outperform both the classic LS and the
state-of-the-art KRF algorithms in terms of NMSE by approximately $10$ dB and
$5$ dB, respectively. The performance gap between the proposed solutions is
small with the ALS algorithm having the best performance in terms of NMSE
only. Regarding the computational complexity, the proposed ALS and HOSVD
solutions have approximately the same complexity as the benchmark ones,
namely, the LS (13) and the KRF [8] methods.

## References

  * [1] N. Rajatheva and et al., “Scoring the terabit/s goal: Broadband connectivity in 6G,” 2020, arXiv:2008.07220.
  * [2] B. Zheng, C. You, W. Mei, and R. Zhang, “A survey on channel estimation and practical passive beamforming design for intelligent reflecting surface aided wireless communications,” IEEE Commun. Surv. Tutor., vol. 24, no. 2, pp. 1035–1071, 2022.
  * [3] S. Gong, X. Lu, D. T. Hoang, D. Niyato, L. Shu, D. I. Kim, and Y.-C. Liang, “Toward smart wireless communications via intelligent reflecting surfaces: A contemporary survey,” IEEE Commun. Surv. Tutor., vol. 22, no. 4, pp. 2283–2314, 2020.
  * [4] C. Pan, G. Zhou, K. Zhi, S. Hong, T. Wu, Y. Pan, H. Ren, M. Di Renzo, A. L. Swindlehurst, R. Zhang, et al., “An overview of signal processing techniques for RIS/IRS-aided wireless systems,” IEEE J. Sel. Topics Signal Processing, vol. 16, no. 5, pp. 883–917, 2022.
  * [5] Y. Lin, S. Jin, M. Matthaiou, and X. You, “Channel estimation and user localization for IRS-assisted MIMO-OFDM systems,” IEEE Transactions on Wireless Communications, vol. 21, no. 4, pp. 2320–2335, 2021.
  * [6] Fazal-E-Asim, A. L. F. de Almeida, B. Sokal, B. Makki, and G. Fodor., “Two-dimensional channel parameter estimation for IRS-assisted networks,” 2023, arXiv:2305.04393.
  * [7] X. Zheng, P. Wang, J. Fang, and H. Li, “Compressed channel estimation for IRS-assisted millimeter wave OFDM systems: A low-rank tensor decomposition-based approach,” IEEE Wireless Communications Letters, vol. 11, no. 6, pp. 1258–1262, 2022.
  * [8] G. T. de Araújo, A. L. de Almeida, and R. Boyer, “Channel estimation for intelligent reflecting surface assisted MIMO systems: A tensor modeling approach,” IEEE J. Sel. Top. Signal Process, vol. 15, no. 3, pp. 789–802, 2021.
  * [9] G. T. de Araújo, P. R. Gomes, A. L. de Almeida, G. Fodor, and B. Makki, “Semi-blind joint channel and symbol estimation in IRS-assisted multi-user MIMO networks,” IEEE Wireless Commun. Lett., 2022.
  * [10] P. R. B. Gomes, G. T. d. Araújo, B. Sokal, A. L. F. d. Almeida, B. Makki, and G. Fodor, “Channel estimation in RIS-assisted MIMO systems operating under imperfections,” IEEE Transactions on Vehicular Technology, pp. 1–14, 2023.
  * [11] K. B. A. Benicio, A. L. F. de Almeida, B. Sokal, Fazal-E-Asim, B. Makki, and G. Fodor., “Tensor-based channel estimation and data-aided tracking in IRS-assisted MIMO systems,” 2023, arXiv:2305.10499.
  * [12] R. W. Heath, N. Gonzalez-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave MIMO systems,” IEEE J. Sel. Top. Signal Process, vol. 10, no. 3, pp. 436–453, 2016.
  * [13] Fazal-E-Asim, F. Antreich, C. C. Cavalcante, A. L. F. de Almeida, and J. A. Nossek, “Two-dimensional channel parameter estimation for millimeter-wave systems using butler matrices,” IEEE Trans. Wireless Commun., vol. 20, no. 4, pp. 2670–2684, 2021.
  * [14] B. Sokal, P. R. Gomes, A. L. de Almeida, and M. Haardt, “Tensor-based receiver for joint channel, data, and phase-noise estimation in MIMO-OFDM systems,” IEEE J. Sel. Top. Signal Process, vol. 15, no. 3, pp. 803–815, 2021.
  * [15] P. Comon, X. Luciani, and A. L. de Almeida, “Tensor decompositions, alternating least squares and other tales,” J. Chemom., vol. 23, no. 7-8, pp. 393–405, 2009.
  * [16] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM review, vol. 51, no. 3, pp. 455–500, 2009.
  * [17] G. Favier, C. A. R. Fernandes, and A. L. de Almeida, “Nested tucker tensor decomposition with application to mimo relay systems using tensor space–time coding (tstc),” Signal Processing, vol. 128, pp. 318–331, 2016.

  *[MIMO]: multiple input multiple output
  *[IRS]: intelligent reflecting surface
  *[B5G]: beyond fifth generation
  *[6G]: sixth generation
  *[SNR]: signal-to-noise ratio
  *[mmWave]: millimeter Wave
  *[BS]: base station
  *[UE]: user equipment
  *[ALS]: alternating least squares
  *[HOSVD]: high order single value decomposition
  *[LS]: least squares
  *[KRF]: Khatri-Rao factorization
  *[AWGN]: additive white Gaussian noise
  *[LOS]: line of sight
  *[NLOS]: non-line of sight
  *[AoA]: angle of arrival
  *[AoD]: angle of departure
  *[SVD]: singular value decomposition
  *[DFT]: discrete Fourier transform
  *[NMSE]: normalized mean square error