#  Terahertz Induced Protein Interactions in a Random Medium

Hadeel Elayan, Andrew W. Eckford, and Raviraj Adve We would like to
acknowledge the support of the National Science and Engineering Research
Council, Canada, through its Discovery Grant program.H. Elayan and R. Adve are
with the Edward S. Rogers Department of Electrical and Computer Engineering,
University of Toronto, Ontario, Canada, M5S 3G4 (e-mail:
<EMAIL_ADDRESS>rsadve@ece.utoronto.ca). A. Eckford is with
the Department of Electrical Engineering and Computer Science, York
University, Ontario, Canada, M3J 1P3 (e-mail: aeckford@yorku.ca). This work
has been presented in part in [1].

###### Abstract

Folding of proteins into their correct native structure is key to their
function. Simultaneously, the intricate interplay between cell movement and
protein conformation highlights the complex nature of cellular processes. In
this work, we demonstrate the impact of Terahertz (THz) signaling on
controlling protein conformational changes in a random medium. Our system of
interest consists of a communication link that involves a nanoantenna
transmitter, a protein receiver, and a channel composed of moving red blood
cells. Due to the system dynamics, we investigate the influence of both the
fast and slow channel variations on protein folding. Specifically, we analyze
the system’s selectivity to asses the effectiveness of the induced THz
interaction in targeting a specific group of proteins under fading conditions.
By optimizing the selectivity metric with respect to the nanoantenna power and
frequency, it is possible to enhance the controllability of protein
interactions. Our probabilistic analysis provides a new perspective regarding
electromagnetically triggered protein molecules, their micro-environment and
their interaction with surrounding particles. It helps elucidate how external
conditions impact the protein folding kinetics and pathways. This results in
not only understanding the mechanisms underlying THz-induced protein
interactions but also engineering these still-emerging tools.

###### Index Terms:

Terahertz signals, protein interactions, probability of folding, fast fading,
slow fading, selectivity, outage.

##  I Introduction

Over the past two decades, the field of nanotechnology has experienced
significant advancements, prompting researchers to explore methods for
establishing reliable communication between nanomachines. This communication
enables nanosensors to independently transmit their sensing data, receive
instructions from a central command center, and collaborate with other
nanomachines when necessary [2]. To achieve nanoscale communication,
scientists have put forward various solutions that consider both molecular and
electromagnetic (EM) communication approaches. From the EM perspective,
plasmonic nano-lasers, plasmonic nanoantennas, and single-photon detectors
have all indicated the significance of the Terahertz (THz) band, which spans
from 0.1 to 10 THz, as a crucial facilitator of communication at the nanoscale
[3].

The viability of utilizing THz intra-body communication has been strengthened
through the characterization of tissues in the THz frequency range [4]. As a
result, scientists have focused on creating models that accurately capture
signal propagation, photo-thermal effects, and sources of noise [5, 6]. For
example, the authors in [7] presented a channel model for THz propagation
inside the human skin. To verify the achieved results, the authors have
compared their findings with THz time-domain spectroscopy measurements of a
skin sample. Moreover, researchers have proposed mechanisms to overcome the
high path-loss and enable wireless communication with deeper implants for
nano-biosensing applications [8]. Recently, the authors in [9] presented an
in-depth overview covering the state of the art in EM nanoscale communication
and networking, where they discussed nanonetworking in the THz band and beyond
while focusing on applications ranging from nano-bio interfaces to quantum
communications.

THz radiation operates within the energy range of molecular interactions,
which encompass hydrogen bonds, intermolecular forces, and the vibration of
macromolecules [10]. This alignment of energy levels elucidates why
biomolecular interactions are highly responsive to emerging THz techniques.
Specifically, proteins, the fundamental components of cells, exhibit
distinctive collective vibrational patterns in the THz frequency range. These
patterns correspond to crucial functional modes, which have been identified as
potential drivers of conformational alterations, ligand binding, and changes
in oxidation state [11].

The dynamics and conformational changes of proteins involve a multitude of
motions driven by various mechanisms, occurring at different timescales and
amplitudes. Protein motions involve a broad range of timeframes, from fast
bond vibrations in the femtosecond range to slow and extensive domain
movements in the millisecond range [12]. These variations have posed
challenges in identifying universal principles governing protein folding [13].
Proteins also experience constant external stresses and crowding from other
biomolecules within cells, membranes, or extracellular spaces [14].

In our previous work, we proposed a hybrid communication model that combines
EM and molecular approaches by utilizing proteins as molecular machines within
the human body [15]. We employed a Markov model to derive the mutual
information and calculate the capacity of the communication system under
different constraints. This analysis was conducted for both a two-state and a
multi-state protein model [16]. In this context, capacity serves as a
quantitative measure that provides researchers with insights into the amount
of information gained by the protein signaling pathway as the protein
transitions from an unfolded to a folded state.

In the aforementioned works, we treated the external nanoantenna force as a
deterministic component, in which the incident EM field directly impacts the
protein population. However, we have not considered the impact of channel
randomness on the EM wave that triggers the desired protein. In a more
realistic model, cells act as obstacles in the path between the transmitter
(Tx) and receiver (Rx). In addition, rapid fluctuations in received signal
strength occur over both short distances and short time periods due to the
mobility of these particles. For instance, in a blood medium, red blood cells
(RBCs) interfere with propagation in a time-varying manner. This phenomenon is
similar to fading in wireless communications, with the impinging signal
suffering from shadowing. The fading can be either slow or fast, depending on
the time constant of the protein dynamics [1].

To bridge the gap, we consider a communication link that consists of a
nanoantenna Tx, a protein Rx, and a channel composed of RBCs. Since a large
number of moving molecules exist between the communicating endpoints, the
channel experiences fading, which randomly attenuates the transmitted
nanoantenna signal. Therefore, we will examine the impact of both fast and
slow fading on our ability to stimulate controlled protein folding behaviors.
To evaluate the system’s performance under fading conditions, we study the
selectivity of the system, which was initially investigated in [17] to
demonstrate the efficiency of the stimulated THz interaction for a protein
population.

The proposed selectivity metric assesses the capability of the nanoantenna to
induce a conformational change within the desired protein population without
affecting untargeted proteins in the network. It requires suppressing the
effect of undesired proteins and enhancing only the impact of targeted ones.
Specifically, in this work, we explore the effect of fading on the ability to
achieve selective protein interactions. The unpredictability of fading
presents a challenge, but by optimizing the selectivity metric with respect to
the nanoantenna power and frequency, it is possible to enhance the
controllability of protein interactions. From this perspective, we make the
following contributions:

  * •

We determine the best fit of the nanoantenna force realizations based on the
developed propagation model.

  * •

We study the phenomenon of fast and slow fading in THz intra-body
communication and its impact on folding dynamics.

  * •

We formulate a joint optimization problem to retrieve the optimal parameters
that maximize the selectivity in fading scenarios.

  * •

We introduce selectivity outage as a metric that indicates the performance of
the slow fading channel.

Studying how the blood channel affects protein dynamics offers a broader
perspective on the functioning and interactions of proteins within a living
organism. This approach helps bridge the gap between in-vitro experiments and
in-vivo conditions. These insights could impact the way experiments are
designed and controlled, as researchers might explore the channel as a tool to
influence folding patterns, investigate specific aspects of folding
mechanisms, or validate hypotheses related to protein folding processes.

The rest of the paper is organized as follows. In Sec. II, we present the
system model for the stimulated protein dynamics. In Sec. III, we present a
propagation model to compute the power received by the protein and calculate
the path-loss experienced due to the shadowing imposed by the particles in the
channel. In Sec. IV, we formulate the expression of the steady-state energy
absorbed by the protein subject to the random nanoantenna force and use it to
modify the Boltzmann distribution. In Sec. V, we introduce the concept of
fading in an intra-body environment and examine its impact on the capability
to invoke controlled interactions in a protein network. In Sec. VI, we
demonstrate our numerical results. Finally, we draw our conclusions in Sec.
VII.

##  II System Model

Our system model involves a nanoantenna Tx and a protein Rx. The channel we
focus on is a blood vessel that contains various intra-body particles. Among
these particles, RBCs are the largest, with a size of approximately 7 microns,
and they are also the most abundant, constituting 45$\%$ of the particles in
human blood. Therefore, RBCs play a crucial role in governing the propagation
of THz waves in the bloodstream. Additionally, the blood cells are surrounded
by blood plasma, which mainly consists of water (92$\%$). We model the plasma
as a lossy medium, characterized by the complex permittivity of water at THz
frequencies [18].

It is noteworthy that in our current work, we consider a dipole nanoantenna
designed for the incident frequency being used, which is tuned to the
vibrational frequency of the desired protein, without delving into the
specific details of the antenna design. For more information on the design of
nanoantennas, we refer interested readers to the works in [3] and [19].

Cells have receptor proteins that bind to signaling molecules and initiate
physiological responses. In our model, we consider our protein as a G-protein-
coupled receptor located on the cell’s exterior. When the receptor protein
receives the nanoantenna signal, it undergoes a conformational change.

Figure 1: System model depicting a communication link: a nanoantenna Tx, a
protein Rx and a channel composed of particles. The model is created with
BioRender.com.

Consider the blood vessel as depicted in Fig. 1. For an $\mathbf{\hat{x}}$
polarized plane wave propagating in the $\mathbf{+\hat{z}}$ direction, the
electric field is given as [20]

$\begin{split}\mathbf{E}(z,t)&=\mathbf{Re}\left[\mathbf{E}(z)e^{-j\omega
t}\right]\mathbf{\hat{x}}\\\ &=E_{o}e^{-\kappa z}\mathbf{Re}\left[{e^{j(\omega
t-\bar{\beta}z)}}\right]\mathbf{\hat{x}}\\\ &=E_{o}e^{-\kappa z}\cos(\omega
t-\bar{\beta}z)\mathbf{\hat{x},}\end{split}$ (1)

where $\mathbf{E}({z})=E_{o}e^{-\varphi z}\mathbf{\hat{x}}$ is the linearly
polarized electric field vector phasor, $\varphi=\kappa+j\bar{\beta}$ is the
complex wave number, $z$ is the propagation distance between the nanoantenna
and the protein, $\kappa$ is the attenuation constant (Nepers/m),
$\bar{\beta}$ is the phase propagation constant (rad/m). Since the human body
is a lossy medium, $\kappa$ and $\bar{\beta}$ have the following general forms
[20]

$\kappa=\omega\sqrt{\left[\frac{\ddot{\mu}\varepsilon^{\prime}_{m}}{2}\left[\sqrt{1+\left(\frac{\varepsilon^{\prime\prime}_{m}}{\varepsilon^{\prime}_{m}}\right)^{2}}-1\right]\right]},$
(2)

$\bar{\beta}=\omega\sqrt{\left[\frac{\ddot{\mu}\varepsilon^{\prime}_{m}}{2}\left[\sqrt{1+\left(\frac{\varepsilon^{\prime\prime}_{m}}{\varepsilon^{\prime}_{m}}\right)^{2}}+1\right]\right]},$
(3)

where $\omega=2\pi f$ is the angular frequency, $\varepsilon^{\prime}_{m}$ and
$\varepsilon^{\prime\prime}_{m}$ are the real and imaginary parts,
respectively, of the complex permittivity of the medium
($\varepsilon^{*}_{m}=\varepsilon^{\prime}_{m}-j\varepsilon^{\prime\prime}_{m}$),
which we consider as a combination of plasma and RBCs, respectively; and
$\ddot{\mu}$ is the relative permeability that acquires a value of one since
most organic materials are non-magnetic at THz frequencies [21].

When a polarized object is subjected to the electric field given in (1), it
induces a dipole moment. If the electric field is inhomogeneous, the field
strength, and therefore the force, acting on each side of the particle will
differ. Dielectrophoresis (DEP) is the force that arises from the interaction
of a dielectric particle with the spatial gradient of an electric field. In
our system, we model the protein as a sphere, which is a common approach for
describing the DEP force acting on a particle [22].

The strength of the induced force, exerted on each end of the dipole, depends
on both the medium and the protein’s electrical properties, on the protein’s
shape and size, as well as on the frequency of the electric field. The time-
averaged force acting on the protein is therefore given as [22]

$f_{o}=2\pi
r^{3}\varepsilon_{m}\mathbf{Re}\left\\{\mathcal{K}\right\\}\nabla|\mathbf{E}_{rms}|^{2},$
(4)

where $\mathbf{E}_{rms}$ is the root mean square value of the electric field
given in (1) and $|\mathbf{E}_{rms}|^{2}=\frac{1}{2}|E^{2}_{o}|e^{-2\kappa
z}$. In addition, $\nabla$ is the mathematical  Del vector operator
($\nabla=\frac{\partial}{\partial x}\hat{x}+\frac{\partial}{\partial
y}\hat{y}+\frac{\partial}{\partial z}\hat{z}$), and
$\nabla|\mathbf{E}_{rms}|^{2}=-\kappa|E^{2}_{o}|e^{-2\kappa z}$. From (4), $r$
is the protein radius, $\varepsilon_{m}=\varepsilon_{o}\varepsilon^{*}_{m}$ is
the absolute permittivity of the medium, and $\mathcal{K}$ is the frequency-
dependent polarization coefficient given by

$\mathcal{K}=\frac{\varepsilon^{*}_{p}-\varepsilon^{*}_{m}}{\varepsilon^{*}_{p}+2\varepsilon^{*}_{m}},$
(5)

where $\varepsilon^{*}_{p}$ is the complex permittivity of the protein sample
[23].

##  III Propagation Model

When an object is illuminated by a wave, the incident power experiences both
absorption and scattering effects caused by particles. These phenomena
contribute to power losses compared to the line-of-sight (LOS) component of
the channel. In our specific propagation scenario, we utilize a first-order
multiple scattering model to describe the propagation of THz waves [24]. This
model accounts for the absorption and scattering that occur along the ray
path. It is particularly applicable when particles primarily exhibit
absorption rather than scattering behavior. This situation arises when the
particle scattering albedo, which represents the ratio of scattering to total
attenuation, is less than 0.5 [25].

###  III-A Signal Transmission

In our scenario, the nanoantenna transmits power $P_{t}$ through a distance
$z$. Since the THz beams are narrow, the beam pattern can be approximated by a
Gaussian function. As such, the nanoantenna beam pattern is given by [26]

$G_{t}(\theta_{t})=\exp\left[-4\ln
2\left(\frac{\theta_{t}}{\theta_{b}}\right)^{2}\right],$ (6)

where $\theta_{t}$ is the angle measured from the horizontal direction, and
$\theta_{b}$ is the 3 dB beamwidth angle. The Rx in our system is a protein.
We treat it as a point Rx with a unity gain, i.e., $G_{r}(\theta_{r})=1$. The
received power is therefore given as

$P_{r}=\frac{P_{t}\lambda_{g}^{2}G_{t}(\theta_{t})G_{r}(\theta_{r})}{(4\pi
z)^{2}}e^{-\gamma}.$ (7)

Here, $\gamma$ is given by

$\gamma=\int_{0}^{z}\rho\left\langle\sigma_{t}\right\rangle\,\,ds,$ (8)

where the density of the molecules $\rho(s)$ and the total cross-section
$\sigma_{t}(s)$ can be functions of the position along the path from the Tx to
the Rx [24]. In specific, $\sigma_{t}=\sigma_{abs}+\sigma_{sca}$ is the sum of
both the absorption and scattering cross-sections, and
$\left\langle\cdot\right\rangle$ represents the average over the size
distribution. The absorption cross-section, $\sigma_{abs}$, and the scattering
cross-section, $\sigma_{sca}$, have been derived in [27].

To take into account the effect of the particle size distribution, we have

$\rho\left\langle\sigma_{t}\right\rangle=\int_{0}^{\infty}\sigma_{t}(x_{d})n(x_{d},\mathbf{r})\,\,dx_{d},$
(9)

where $n(x_{d},\mathbf{r})$ is the number of RBCs per unit volume located at
$\mathbf{r}$ having a range of sizes, in this case, the diameter of the RBC,
between $x_{d}$ and $x_{d}+dx_{d}$. Therefore, the number density distribution
as a function of the particle size follows a log-normal distribution given by
[28]

$n(x_{d},\mathbf{r})=\frac{1}{x_{d}\sqrt{(2\pi)}\ln(\sigma_{g})}e^{-\frac{(\ln(x_{d})-\ln(r_{g}))^{2}}{2\ln(\sigma_{g})^{2}}},$
(10)

where $r_{g}$ is the mean geometric diameter of the scatterers, and
$\sigma_{g}$ is the corresponding standard deviation. We note that skewed
distributions are particularly common when values cannot be negative due to
the physical aspects of the problem. In (7),
$\lambda_{g}=\lambda/n^{\prime}(\omega)$, where $n^{\prime}(\omega)$ is the
real part of the protein refractive index. We compute $n^{\prime}(\omega)$
using the protein permittivity expression given in [23].

###  III-B Path-Loss and Shadowing

The dynamic nature of intra-body communications results in random shadowing
events. These events occur when the transmitted EM waves encounter obstacles,
in this case, particles, that obscure the direct LOS signal path. The severity
of shadowing depends on factors such as the size and material properties of
the obstructing objects, the frequency of the signal, and the distance between
the Tx and Rx.

According to the central limit theorem, the attenuation, $X_{\sigma}$, can be
modeled as a Gaussian random variable in the log-scale, referred to as log-
normal shadowing. Hence, $X_{\sigma}$ has the following distribution

$f(X_{\sigma})=\frac{1}{\sqrt{2\pi\sigma_{X}^{2}}}e^{-\frac{(X_{\sigma}-\mu_{X})^{2}}{2\sigma_{X}^{2}}},$
(11)

where $\mu_{X}$ and $\sigma^{2}_{X}$ are the mean and variance in dB of the
shadowing variable, respectively. In this paper, the parameters defining the
pdf of $X_{\sigma}$ are extracted by implementing Monte Carlo simulations
[27]. By incorporating the effect of shadowing and knowing the power received
from (7), the path-loss can be formulated as

$\left[L_{LOS}\right]_{\mathbf{dB}}=\left[\frac{(4\pi
z)^{2}}{\lambda^{2}_{g}G_{t}(\theta_{t})G_{r}(\theta_{r})}e^{\gamma}\right]_{\mathbf{d}\mathbf{B}}+X_{\sigma}.$
(12)

##  IV Modeling Protein Dynamics

###  IV-A Langevin Equation

To model protein dynamics, we analyze the forces acting on the system. These
forces include impacts from the external force exerted by the nanoantenna,
which is subject to channel randomness, as well as forces arising from the
Brownian motion of the surrounding particles. By solving the Langevin
equation, researchers can gain insights into the behavior and dynamics of
proteins, including their diffusion properties, conformational changes, and
interactions with their surroundings. Therefore, we deploy the Langevin
stochastic equation, driven by these forces as follows [29]

$m\frac{d^{2}x}{dt^{2}}+\ddot{\beta}\frac{dx}{dt}+kx(t)=f_{ex}(t)+f_{\zeta}(t),$
(13)

where $x=x(t)$ is the protein coordinate, $m$ is the protein mass, $k$ is the
protein spring constant, and $\ddot{\beta}$ is the damping coefficient. In
addition, $f_{ex}(t)$ is the external driving nanoantenna force, which can be
represented as

$f_{ex}(t)=f_{o}\cos(\omega_{a}t+\Theta),$ (14)

while $f_{\zeta}(t)$ is the internal force acting on the protein due to the
random motion of the particles suspended in the medium.

####  IV-A1 Nanoantenna External Force

To relate the nanoantenna external force, $f_{ex}(t)$, to the propagation
model, we first exploit the relationship between the received power, $P_{r}$,
absorbed by the protein and the impinging nanoantenna electric field strength,
$E_{o}$. This is given as [20]

$P_{r}=\frac{1}{2\mathbf{R}\mathbf{e}\left\\{\eta^{*}_{c}\right\\}}E_{o}^{2}A_{r},$
(15)

where $\eta^{*}_{c}$ is the intrinsic impedance in a lossy medium, given by
[20]

$\eta_{c}=\sqrt{\frac{\ddot{\mu}}{\varepsilon^{\prime}_{m}}}\left(1-j\frac{\varepsilon^{\prime\prime}_{m}}{\varepsilon^{\prime}_{m}}\right)^{-\frac{1}{2}}.$
(16)

In addition, $A_{r}$ is the receiving cross-section for the THz wave incident
on the protein. It is given as

$A_{r}=\frac{\lambda_{g}^{2}}{4\pi}G_{r}(\theta_{r}).$ (17)

We then find the relationship between the incident electric field and the
force stimulating the protein by solving for $E^{2}_{o}$. This is done by
equating (7) to (15), which yields

$E^{2}_{o}=\frac{2P_{t}G_{t}(\theta_{t})G_{r}(\theta_{r})\mathbf{R}\mathbf{e}\left\\{\eta^{*}_{c}\right\\}}{L_{LOS}\,A_{r}}.$
(18)

Finally, we substitute (18) in (4) to obtain

$f^{2}_{o}=\left(\frac{-4\pi
r^{3}\varepsilon_{m}\mathbf{Re}\left\\{\mathcal{K}\eta^{*}_{c}\right\\}\kappa
P_{t}G_{t}(\theta_{t})G_{r}(\theta_{r})e^{-2\kappa
z}}{L_{LOS}A_{r}}\right)^{2}.$ (19)

We note here that the randomness in the nanoantenna force stems from the path-
loss, $L_{LOS}$, experienced by the channel due to the large number of
molecules.

Using COMSOL Multiphysics® software, we built a physical model that mimics our
intra-body scenario and fitted $f^{2}_{o}$ as demonstrated in Fig. 2. Among
the candidate distributions, the gamma distribution was the best fit, as
indicated by the log-likelihood value provided in Table I. The log likelihood
criterion is a measure of goodness of fit for any model, in which the higher
the log-likelihood value, the better the model.

Figure 2: Fitting $f^{2}_{o}$ values obtained from the received EM field strength in COMSOL, where $\alpha=2.5\times 10^{-17}$ and $\beta=1.9$. TABLE I: Comparison of distributions for the best-fit of $f^{2}_{o}$ Distribution |  Best-Fit |  Log likelihood
---|---|---
Gamma |  1 |  7.167e+04
Lognormal |  2 |  7.163e+04
Weibull |  3 |  7.161e+04
Nakagami |  4 |  7.155e+04
Exponential |  5 |  7.143e+04
Rayleigh |  6 |  7.135e+04

Thereby, $f^{2}_{o}$ follows a gamma distribution given by

$g(f_{o}^{2};\alpha,\beta)=\frac{f_{o}^{2(\alpha-1)}e^{-\beta
f^{2}_{o}}\beta^{\alpha}}{\Gamma(\alpha)}.$ (20)

The gamma distribution is a two-parameter family of continuous probability
distributions, comprising a shape parameter $\alpha$ and a rate parameter
$\beta$. Changing $\alpha$ alters the shape of the probability density
function, while adjusting $\beta$ either stretches or compresses the range of
the gamma distribution. It’s worth noting that varying the density, $\rho$,
and size, $x_{d}$, of the RBCs results in different values for $\alpha$ and
$\beta$.

When $\alpha=1$, the resulting gamma distribution is equivalent to the
exponential distribution. Interestingly, single exponential decays serve as a
natural starting point for analyzing protein relaxation phenomena, as they are
phenomenological descriptors of simple dashpots—dampers that exhibit
frictional resistance to motion. In addition, the distribution of protein
abundances in a cell is well described in the literature by a gamma
distribution [30].

Moreover, $\Theta$ in (14) is a random variable uniformly distributed between
$\left[-\pi,\pi\right]$, and $\omega_{a}$ is the nanoantenna frequency. Since
the amplitude and phase in (14) are independent, we can prove that $f_{ex}(t)$
is a wide-sense stationary process, by computing its mean and autocorrelation
function. For a gamma distribution, the mean and variance are given by
$E[f^{2}_{o}]=\frac{\alpha}{\beta}$ and
$Var[f^{2}_{o}]=\frac{\alpha}{\beta^{2}}$, respectively. The auto-correlation
function, for time-lag $\tau$, yields

$R_{f_{ex}}(\tau)=\frac{\alpha}{2\beta}\cos(\omega_{a}\tau).$ (21)

####  IV-A2 Force due to Brownian Motion

In (13), the internal stochastic force $f_{\zeta}(t)$ is governed by a white-
noise fluctuation-dissipation relation as follows [31]

$f_{\zeta}(t)={\sqrt{\Gamma}}\zeta(t).$ (22)

Here, $\zeta(t)$ is a white, Gaussian, random process with moments [31]

$\begin{split}\langle\zeta(t_{1})\rangle&=0\\\
\langle\zeta(t_{1})\zeta(t_{2})\rangle&=\delta(t_{1}-t_{2}).\end{split}$ (23)

In (22), $\Gamma=2\bar{\gamma}mk_{b}T$ ($k_{b}$ is Boltzmann’s constant and
$T$ is the temperature in Kelvin) denotes the strength of the noise
$\zeta(t)$. It fixes the amplitude of the fluctuation in the random force in
terms of both the temperature and the dissipation coefficient [31].

The Green’s function is a powerful mathematical tool to solve inhomogeneous
differential equations. For an arbitrary forcing term, the solution to the
equation is formulated by integrating the Green’s function against the forcing
term. To solve (13), we first re-write it as

$\frac{d^{2}x}{dt^{2}}+\bar{\gamma}\frac{dx}{dt}+\omega_{o}^{2}x(t)=\frac{1}{m}\left[f_{ex}(t)+f_{\zeta}(t)\right],$
(24)

where $\bar{\gamma}={\ddot{\beta}}/{m}$ and $\omega_{o}^{2}=k/m$ [32]. Here,
$\omega_{o}$ is the natural frequency of the protein, and $\bar{\gamma}$ is
the damping constant that governs the magnitude of the protein vibrational
resonances.

In the context of modeling proteins, the interatomic forces are often
approximated using harmonic potentials. This approximation is based on the
assumption that the potential energy between atoms within the protein can be
described as a harmonic oscillator potential. In fact, in [33], the author
showed that the vibrational modes of protein molecules are not significantly
modified when interactions are replaced by Hookean springs, for all atom pairs
whose distance is smaller than a given cut-off value.

The frequency-domain solution of the Green’s function corresponding to (24) is
found as [34]

$G(\omega)=\frac{1}{m\left[-\omega^{2}-j\bar{\gamma}\omega+\omega^{2}_{o}\right]},$
(25)

where $G(\omega)$ represents the transfer function of the protein dynamics.
The solution of (24) can be decomposed into two parts, one related to the
nanoantenna external force $f_{ex}(t)$ and the other is related to the
stochastic force $f_{\zeta}(t)$, resulting in
$x(t)=~{}x_{ex}(t)+x_{\zeta}(t)$.

In our model, we consider $\omega_{o}>\frac{\bar{\gamma}}{2}$. Physically,
this condition satisfies protein collective vibrations, where protein modes in
the THz frequency range have been shown to be underdamped even in aqueous
solutions [35, 36].

###  IV-B Energy of Driven Protein

Proteins are dynamic entities capable of adopting various conformations to
perform their biological functions. These conformational changes are often
triggered by the exchange of energy between the protein and its environment.
The requirement for energy transfer to induce a conformational change
highlights the importance of connecting THz band signals with protein
molecules to initiate resonance. When an external harmonic excitation has a
frequency that matches one of the natural frequencies of the system, resonance
occurs, leading to an increase in the vibrational amplitude of the structure
[32].

In protein folding, the energy landscape and the resulting folding pathway
play a crucial role in determining the probability of a protein adopting its
native conformation. To assess the impact of THz signaling on protein folding
probability, we must first calculate the total energy absorbed by the protein
when subjected to both the nanoantenna force and the force resulting from the
surrounding Brownian motion. To do so, we sum the kinetic and potential energy
contents [37]

$\left\langle E_{tot}\right\rangle_{ss}=\frac{1}{2}m\left\langle
v^{2}(t)\right\rangle_{ss}+\frac{1}{2}k\left\langle
x^{2}(t)\right\rangle_{ss}.$ (26)

We denote by $\left\langle\cdot\right\rangle_{ss}$ the statistical average
value in steady-state, where all transient effects die out.

To evaluate the potential energy, we find the average steady-state squared
displacement of the protein, $\left\langle x^{2}(t)\right\rangle_{ss}$, which
is decomposed into

$\left\langle x^{2}(t)\right\rangle_{ss}=\left\langle
x_{ex}^{2}(t)\right\rangle_{ss}+\left\langle
x_{\zeta}^{2}(t)\right\rangle_{ss}.$ (27)

$\left\langle x^{2}(t)\right\rangle_{ss}$ is a measure of the spatial extent
of the protein motion. As detailed in [17], $\left\langle
x_{\zeta}^{2}(t)\right\rangle_{ss}$ yields $\frac{k_{b}T}{k}$. In terms of
$\left\langle x_{ex}^{2}(t)\right\rangle_{ss}$, we start by examining
$f_{ex}(t)$ given in (14). Since the latter is a wide-sense stationary random
process with autocorrelation function $R_{f_{ex}}(\tau)$, its power spectral
density (PSD) $S_{{ex}}(\omega)$ is introduced as the Fourier transform of
$R_{f_{ex}}(\tau)$. This relationship is given by

$\begin{split}S_{{ex}}(\omega)&=\int_{-\infty}^{\infty}R_{f_{ex}}(\tau)e^{-j\omega\tau}d\tau\\\
&=\frac{\pi\alpha}{2\beta}\left[\delta(\omega-\omega_{a})+\delta(\omega+\omega_{a})\right].\end{split}$
(28)

Using the PSD of the driving force and the frequency-domain solution of
Green’s function in (25), the output PSD is given by

$S_{x}(\omega)=|G(\omega)|^{2}\cdot S_{{ex}}(\omega).$ (29)

Finally, from (29), the average squared displacement as a result of the
nanoantenna force, $\left\langle x^{2}_{ex}(t)\right\rangle_{ss}$, is found as

$\begin{split}\left\langle
x^{2}_{ex}(t)\right\rangle_{ss}&=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{x}(\omega)d\omega\\\
&=\frac{1}{2\pi}\frac{\pi\alpha}{2\beta}\int_{-\infty}^{\infty}\frac{(\delta(\omega-\omega_{a})+\delta(\omega+\omega_{a}))}{m^{2}[(\omega_{o}^{2}-\omega^{2})^{2}+(\omega\bar{\gamma})^{2}]}d\omega\\\
&=\frac{\alpha}{2\beta
m^{2}[(\omega_{o}^{2}-\omega_{a}^{2})^{2}+(\omega_{a}\bar{\gamma})^{2}]}.\end{split}$
(30)

To calculate the kinetic energy, we need to compute the average steady-state
squared velocity of the protein, $\left\langle v^{2}(t)\right\rangle_{ss}$,
which is decomposed into

$\left\langle v^{2}(t)\right\rangle_{ss}=\left\langle
v_{ex}^{2}(t)\right\rangle_{ss}+\left\langle
v_{\zeta}^{2}(t)\right\rangle_{ss}.$ (31)

From the derivation in [17], we know that $\left\langle
v_{\zeta}^{2}(t)\right\rangle_{ss}$ yields $\frac{k_{b}T}{m}$. Further, the
steady-state average squared velocity resulting from the nanoantenna force can
be found from the Fourier transform of the derivative of the displacement.
From (30), we have

$\begin{split}\left\langle
v^{2}_{ex}(t)\right\rangle_{ss}&=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{x}(\omega)\omega^{2}d\omega\\\
&=\frac{1}{2\pi}\frac{\alpha\pi}{2\beta}\int_{-\infty}^{\infty}\frac{\omega^{2}(\delta(\omega-\omega_{a})+\delta(\omega+\omega_{a}))}{m^{2}[\left(\omega_{o}^{2}-\omega^{2})^{2}+(\omega\bar{\gamma}\right)^{2}]}d\omega\\\
&=\frac{\alpha\omega^{2}_{a}}{2\beta
m^{2}[(\omega_{o}^{2}-\omega_{a}^{2})^{2}+(\omega_{a}\bar{\gamma})^{2}]}.\end{split}$
(32)

Finally, we find the total steady-state energy of the driven protein motion by
substituting (27) and (31) in (26) to yield

$\begin{split}\left\langle
E_{tot}\right\rangle_{ss}=&\underbrace{\frac{\alpha}{4\beta
m[(\omega_{o}^{2}-\omega_{a}^{2})^{2}+(\omega_{a}\bar{\gamma})^{2}]}(\omega_{a}^{2}+\omega^{2}_{o})}_{\text{THz
Contribution}}+\underbrace{k_{b}T}_{\text{Noise Effect}}.\end{split}$ (33)

If we compare $\left\langle E_{\text{tot}}\right\rangle_{\text{ss}}$ in (33)
to our previous work in [17], we can observe the influence of the gamma
distribution parameters, specifically $\alpha$ and $\beta$, on the energy
expression. $\alpha$ represents the amplitude of the squared force, while
$\beta$ governs the extent of path-loss. This clearly underscores the
influence of the channel on the incident nanoantenna force, with the values of
$\alpha$ and $\beta$ determining the amount of energy absorbed by the protein
particle.

###  IV-C Boltzmann Distribution

Energy functions are linked to molecular conformation through the Boltzmann
distribution. This distribution provides the probability of a system being in
a specific state based on the energy of that state and the system’s
temperature [38]. In the context of protein folding, the Boltzmann
distribution aids in understanding the equilibrium populations of different
folding states. At a given temperature, protein states with lower energies are
more likely to be populated, while states with higher energies have lower
probabilities. The state with the lowest energy corresponds to the protein’s
native conformation, while higher energy states represent various unfolded or
partially folded conformations.

In this work, we model the protein as having two states, unfolded and folded.
Through the folding structure, the protein acquires a conformation that is
biologically functional. The rate of protein folding is given by [39]

$r_{f}=r_{0}\exp\left(\frac{-E_{f}}{k_{b}T}\right),$ (34)

where $E_{f}$ denotes the free-energy associated with the folded protein
state, and $r_{0}$ is a scale factor which preserves detailed balance. In our
case, we amend (34) to incorporate (33) as follows

$r_{f}=r_{0}\exp\left(\frac{-E_{f}+\langle
E_{tot}\rangle_{ss}}{k_{b}T}\right).$ (35)

The unfolding transition is a relaxation process, which returns the protein to
the unfolded state. This process is considered independent of the imposed
signal [40] and is given by

$r_{u}=r_{0}\exp\left(\frac{-E_{u}}{k_{b}T}\right),$ (36)

where $E_{u}$ denotes the free-energy associated with the unfolded protein
state. Consequently, the rate of change in the protein folding probability is
given by

$\frac{d}{dt}p_{f}(t)=-r_{u}p_{f}(t)+r_{f}(1-p_{f}(t)).$ (37)

Here, $p_{f}(t)$ is the probability that the protein is in the folded state at
time $t$ (the corresponding probability of the unfolded state is given by
$p_{u}(t)=1-p_{f}(t)$). From (37), the steady-state solution can be found as

$p_{F}=\frac{r_{f}}{r_{f}+r_{u}}.$ (38)

Substituting (35) and (36) in (38) yields

$p_{F}=\frac{1}{1+\exp\left(\frac{\Delta E-\langle
E_{tot}\rangle_{ss}}{k_{b}T}\right)},$ (39)

where $\Delta E=E_{f}-E_{u}$ is the protein free-energy.

The state of a single protein can be regarded as a Bernoulli random variable
with a probability of success, $p_{F}$. When considering a system composed of
$n$ proteins, the random number of folded proteins, $n_{F}$, follows a
binomial distribution. Due to the large diversity of proteins seen in an
intra-body environment, the binomial distribution can be approximated by a
normal distribution with mean, $\mu=np_{F}$, and variance,
$\sigma^{2}=np_{F}(1-p_{F})$. Finally, we define protein population as the
number of copies of a protein molecule in a cell. Hence, the number of folded
proteins of each protein population, $i$, can be expressed as a normal
distribution as follows

$n_{F,i}\sim\mathcal{N}\left(np_{F,i},np_{F,i}(1-p_{F,i})\right).$ (40)

##  V The Impact of Fading on THz-Induced Protein Interactions

Fading refers to the attenuation of the THz nanoantenna signal caused by the
presence of RBCs. Furthermore, due to the mobility of RBCs, fading varies over
time. This variation has implications for controlling the dynamics of the
desired protein population and places certain requirements on the system
parameters. To understand the impact of the channel on the receiver, we need
to determine how quickly the communication channel changes relative to the
protein folding time. Therefore, we use the coherence time, denoted as
$T_{c}$, as a metric to characterize the nature of the channel in the time
domain [41]. $T_{c}$ is inversely proportional to the maximum Doppler spread,
$f_{m}$. We have $T_{c}\approx\frac{1}{f_{m}}$, where $f_{m}$ is computed
using [41]

$f_{m}=\frac{v}{c}f_{c}.$ (41)

Here, $f_{c}$ is the center frequency of the nanoantenna (in our case it is
the resonant frequency of the protein of interest), $v$ is the velocity of
RBCs, and $c$ is the speed of light.

Fast fading occurs when the channel impulse response changes rapidly within
the symbol duration $T_{s}$. In this case, high Doppler spread is observed,
where $T_{c}\ll T_{s}$ [41], and $T_{s}$ represents the protein folding time.
Since fast fading is characterized by the rapid fluctuations of the signal
over small distances, it is considered an ergodic process that allows the use
of average metrics for its assessment. Therefore, to compute the probability
of the protein in a folded state, we use the form of the probability $p_{F}$
given in (39). Basically, $\langle E_{tot}\rangle_{ss}$ in (39) represents the
averaged energy, which captures the incident wave randomness that is changing
very fast.

On the other hand, slow fading arises when the channel impulse response is
considered roughly constant over the period of use [41]. In slow fading, low
Doppler spread is observed, where $T_{c}\gg T_{s}$. Slow fading is a long-term
fading effect that changes the mean value of the received signal. To compute
the probability of the protein in a folded state under such conditions, we
take the expectation of $p_{F}$ as

$\begin{split}E[p_{F}]&=\int_{0}^{\infty}\frac{1}{1+\exp(\frac{\Delta
E-E_{i}}{k_{b}T})}g(E_{i};\alpha,\beta)\,\,dE_{i},\end{split}$ (42)

where the instantaneous energy $E_{i}$ follows a gamma distribution,
$g(E_{i};\alpha,\beta)$, stemming from the fact that $E_{i}\propto f^{2}_{o}$.

In (42), $E_{i}=f_{o}\times\Delta x$ results from the interaction between the
nanoantenna force, $f_{o}$, and the protein molecule, where $\Delta x$
corresponds to a conformational change in the protein structure. From the
relationship between $f_{o}$, $\Delta x$ and the protein stiffness $k$, we can
write $\Delta x=f_{o}/k$. Finally, from the vibrational frequency of the
protein structure, we have $k=m\omega^{2}_{o}$ [32]. The instantaneous energy
therefore yields

$E_{i}=\frac{f^{2}_{o}}{m\omega^{2}_{o}},$ (43)

where $f^{2}_{o}$ is given in (19).

###  V-A Selectivity in Fading Scenarios

Selectivity refers to the nanoantenna’s ability to induce a conformational
change in the desired protein population without affecting the conformation of
other untargeted proteins in the system. It is a metric that assesses the
controllability of THz signals over protein networks. In our previous work
[17], we proposed the selectivity metric as

$S({\mu_{d}},\mu_{ud})=\frac{\mu_{d}-\mu_{ud}}{\max(\mu_{d},\mu_{ud})},$ (44)

where $\mu_{d}=n_{d}p_{F,d}$ and $\mu_{ud}=n_{un}p_{F,ud}$ are the means of
the desired and undesired protein populations, respectively. The developed
selectivity metric has several properties which makes it a powerful tool for
evaluating protein interactions in the system as discussed in [17]. Most
importantly, as the metric relies on the stationary probability difference, it
encompasses both the system’s mechanical parameters and the energy stored in
the THz-stimulated protein.

Nevertheless, the presence of fading has an impact on the system’s selectivity
by potentially causing the reception of a weak signal. This weak signal may
lead to a temporary disruption in the communication link between the
nanoantenna and the target protein due to a significant decrease in power.
This, in turn, affects the system’s ability to target the desired protein
population. Although fading in an intra-body scenario is characterized by its
unpredictability and stochastic nature, our metric can still be optimized to
ensure that the system achieves the highest level of selectivity possible.
Understanding the interactions between proteins and components within the
bloodstream, and how these interactions influence disease progression, can aid
in identifying potential therapeutic targets. By concentrating on proteins
involved in specific pathways, researchers can develop therapies that
selectively modulate processes associated with the disease.

In our model, the nanoantenna Tx power, $P_{t}$, and frequency, $\omega_{a}$
are the only system parameters that can be controlled. We, therefore,
formulate a joint optimization problem to maximize the selectivity with
respect to $P_{t}$ and $\omega_{a}$. In a fast-fading scenario, the
optimization problem yields

$\displaystyle\underset{P_{t},\text{$\omega_{a}$}}{\text{max.}\,\,}$
$\displaystyle\text{$S(\mu_{d},\mu_{ud})$}$ (45) subject to $\displaystyle
0\leq p_{F,d}\leq\ 1$ $\displaystyle 0\leq p_{F,ud}\leq\ 1.$

The optimization problem is first solved by expressing $p_{F}$ in terms of
$P_{t}$ and $\omega_{a}$. Then, using the fixed point iteration method [42],
we numerically solve the joint optimization problem in (45), where the optimal
values that maximize the selectivity are retrieved.

###  V-B Probability of Selectivity Outage

In the case of slow fading, the selectivity is given as

$\begin{split}S(E[\mu_{d}],E[\mu_{ud}])&=\frac{E[\mu_{d}]-E[\mu_{ud}]}{\max(E[\mu_{d}],E[\mu_{ud}])},\end{split}$
(46)

where $E[\mu]=nE[p_{F}]$ and $E[p_{F}]$ is solved using (42). The optimization
problem is formulated as

$\displaystyle\underset{P_{t},\text{$\omega_{a}$}}{\text{max.}\,\,}$
$\displaystyle\text{$S(E[\mu_{d}],E[\mu_{ud}])$}$ (47) subject to
$\displaystyle 0\leq p_{F,d}\leq\ 1$ $\displaystyle 0\leq p_{F,ud}\leq\ 1.$

Similar to the previous case, (47) is numerically solved using fixed point
iteration and the optimal $P_{t}$ and $\omega_{a}$ are found.

It’s important to note that in conditions of slow fading, the wireless channel
is categorized as non-ergodic, meaning that the capacity of the channel
becomes a random variable [41]. In such scenarios, a measure of the channel’s
quality is the  outage probability, which represents the likelihood of the
information rate falling below the required threshold rate. To maintain
acceptable communication performance, a minimum signal level is necessary,
ensuring that the received signal maintains sufficient strength during “non-
fade intervals”. However, if the signal level drops below this threshold, the
system will experience insufficient signal strength, resulting in what is
commonly referred to as an “outage.”

The concept of ergodicity is often employed in the context of protein dynamics
to describe the statistical properties of conformational changes over time.
However, certain protein systems exhibit non-ergodic behavior, indicating that
their conformational dynamics cannot be easily sampled or described by a
single trajectory or ensemble of conformations. Non-ergodicity and slow fading
are related in the sense that slow fading can lead to non-ergodic behavior in
wireless communication channels. To study slow fading in an intra-body
scenario, we define selectivity outage as a situation where the selectivity
falls below a certain threshold, i.e., $S(E[\mu_{d}],E[\mu_{ud}])<\gamma_{o}$.
Selectivity can drop below the defined threshold when either the desired
protein population receives a minimal amount of power, insufficient for
folding the targeted proteins, or when the system receives a substantial
amount of power, causing all proteins to fold. In both cases, we lose the
ability to target a specific protein.

In this context, we introduce the probability of selectivity outage as the
probability of attaining a selectivity below the defined threshold, i.e.,
$\mathrm{Pr}(S(E[\mu_{d}],E[\mu_{ud}])<\gamma_{o})$. A trade-off must exist
between the selectivity threshold and the outage probability since the higher
the threshold, the higher the probability of selectivity outage. When an
outage occurs in the system, the nanoantenna must re-transmit the signal. This
not only results in unconstrained delays but also sacrifices the system’s
energy by consuming more power. Different medical applications often require
varying degrees of selectivity depending on the specific goals and
requirements of the application. Therefore, the application must dictate
whether or not the system can sustain being in an outage.

The probability of selectivity outage serves as a metric for assessing the
quality of nanoantenna transmission in an intra-body network under slow fading
conditions. By setting a threshold value, the best performance achievable by
the nanoantenna is to maintain a selectivity level equal to that threshold.
Reliable controllability is achieved when this threshold is met; otherwise, an
outage occurs.

##  VI  Results

TABLE II: Simulation Parameters Parameters |  Value |  Ref.
---|---|---
Rhodopsin resonant frequency |  1.36 THz |  [43]
Rhodopsin stiffness |  3 N/m |  [44]
Rhodopsin mass |  $1.62\times 10^{-24}$ kg |  -
Bacteriorhodopsin resonant frequency |  1.13 THz |  [43]
Bacteriorhodopsin stiffness |  1.9 N/m |  [44]
Bacteriorhodopsin mass |  $1.48\times 10^{-24}$ kg |  -
Damping constant ($\bar{\gamma})$ |  0.3 THz |  [35]
Free energy ($\Delta E$) |  6 kbT |  [45]
Tx Power ($P_{t}$) |  35 $\mu$W |  -
Tx Gain ($G_{t}$) |  1.5 (1.76 dBi) |  -
Protein radius ($r$) |  5 nm |  -
Particle density ($\rho$) |  1000 particles/mm3 |  -

###  VI-A COMSOL Model

We constructed a physical system that replicates our intra-body scenario using
COMSOL Multiphysics®. For our analysis, we utilized the simulation parameters
provided in Table II. Specifically, our model focuses on the transmission of
an EM wave from a nanoantenna to the protein population of interest within a
segment of a blood vessel. In our representation, we employed the
Electromagnetic Waves, Frequency domain interface to depict the nanoantenna as
a short dipole, which is a realistic assumption for short-range communication
between nanomachines [18]. We configured the nanoantenna’s frequency to match
the vibrational frequency of the target protein, rhodopsin (i.e., $1.36$ THz).
The protein’s radius was set at $5$ nm, and we considered a population of 1000
proteins arranged in a cylindrical shape. Furthermore, we incorporated random
particles that mimic RBCs into the field. These particles have a radius of
$4~{}\mu$m and dielectric properties as described in [46].

Fig. 3 illustrates the electric field strength received by the targeted
protein population at a distance of $z=500$ $\mu$m, ensuring that we are
operating in the far-field. From the obtained electric field values, we
calculated the force impinging on the protein particles using (19). Then, by
sampling the force realizations every $2$ ms, we found the distribution of the
force that best fits the data based on a maximum likelihood estimation
approach.

Figure 3: Electric field pattern in V/m for propagation at 1.36 THz.

The objective of this work is to demonstrate the influence of a channel
composed of a large number of moving particles on protein folding probability.
This is accomplished by determining the values of both the steady-state
energy, $\langle E_{tot}\rangle_{ss}$, and the instantaneous energy, $E_{i}$,
required to compute $p_{F}$ in (39) and $E[p_{F}]$ in (42), respectively. To
create a dynamic environment, we utilized the Particle Tracing for Fluid Flow
module in COMSOL. We varied the velocity of RBCs to establish a time-varying
channel within the range of $0$-$150$ milliseconds (ms). Notably, according to
the Protein Folding Kinetics Database, the folding timescales of most proteins
are on the order of milliseconds, with a median of $5$ ms for two-state
proteins [47]. Therefore, we set the protein folding time in our system to
$T_{s}=5$ ms. Then, we calculated $T_{c}$ as summarized in Table III. It’s
worth noting that the speed of RBCs reaches approximately $2$-$3$ mm/s in
brain capillaries, $10$ mm/s in brain arterioles, and exceeds $20$-$30$ mm/s
in larger pial vessels. Standard scanning systems typically provide scanning
velocities ranging from $5$ to $20$ mm/s [48].

TABLE III: Fading scenarios in an intra-body environment. Scenario |  Velocity of RBCs |  Channel Coherence Time
---|---|---
Fast Fading |  200 mm/s |  $T_{c}=1$ ms
Slow Fading |  5 mm/s |  $T_{c}=44$ ms

Fig. 4 shows the received power in dB versus time for a single EM signal
realization under fast and slow fading scenarios. These scenarios are achieved
by varying the velocity of RBCs. In fact, every time the nanoantenna transmits
an EM signal, the strength of the impinging electric field is computed, and
accordingly, the received power is calculated using (15).

Figure 4: Power received in dB by the protein under fast and slow fading
conditions.

In addition, Figs. 5(a) and 5(b) illustrate the mean number of folded proteins
($\mu=np_{F}$) versus time for the fast and slow fading scenarios, where
$n=1000$. Under fast fading, we can notice from Fig. 5(a) the impact of the
fluctuating channel response on the protein folding probability. Since fast
fading is an ergodic process, we are interested in the time-average value,
represented by the orange line in Fig. 5(a). Consequently, from Fig. 5(a), the
average number of proteins in the folded state is $\mu=np_{F}=772$.

In Fig. 5(b), we observe the effect of slow fading, characterized by a roughly
constant amplitude and phase on the number of proteins in the folded state. In
this case, the number of folded proteins averaged over time is
$\mu=nE[p_{F}]=848$. The lower number of folded proteins observed under fast
fading conditions, compared to the slow fading case, can be attributed
biologically to the velocity of RBCs. The high velocity limits the cellular
volume accessible to a protein polypeptide chain, thereby affecting its
ability to attain a folded state. While one might argue that the velocity of
RBCs in the fast fading example is high compared to real-world scenarios,
different combinations of RBC and protein timescales can lead to fast-fading
scenarios. Therefore, the provided example offers a clear visualization of its
impact.

While the mean number of folded proteins is higher in the slow fading case,
there is a noticeable drop in the number of folded proteins between $110$ and
$150$ ms. This decline can be attributed to the fact that slow fading may
affect the protein’s configuration, causing it to become trapped in an
intermediate state within a restricted region of phase space. In other words,
the protein may become stuck in an intermediate state, preventing it from
achieving a fully folded conformation. In such a scenario, the system would
require the nanoantenna to re-transmit the signal, resulting in uncontrolled
delays and stringent power requirements on the link between the nanoantenna
and the targeted protein population.

Figure 5: Mean number of proteins in a folded state (a) under fast fading
scenario. (b) under slow fading scenario.

###  VI-B Selectivity under Fading Scenarios

To further investigate this phenomenon, we examine the impact of slow fading
on the system’s selectivity. Therefore, we introduce another protein
population, bacteriorhodopsin, in the vicinity of the targeted rhodopsin
population to simulate a realistic protein network. The vibrational frequency
of bacteriorhodopsin and its mechanical properties can be found in Table II.
To calculate the selectivity, we use (44), where we initially compare the mean
number of folded proteins under slow fading conditions for both populations
within the network, as illustrated in Fig. 6(a). It’s worth noting that the
nanoantenna frequency is tuned to the vibrational frequency of the rhodopsin
(targeted) protein population.

Fig. 6(b) demonstrates the selectivity of protein rhodopsin, with the orange
line representing the threshold we have set for this example,
$\gamma_{o}=0.6$. As expected, we can clearly observe a correlation between
the performance in Fig. 6(a) and Fig. 6(b). In cases where the mean number of
folded proteins of both populations is comparable, the selectivity is low,
while in cases where the the mean number of folded proteins in both
populations differs significantly, the selectivity is high.

Figure 6: (a) Mean number of folded proteins of rhodopsin (targeted) and
bacteriorhodopsin. (b) Selectivity of rhodopsin ($\gamma_{o}$= 0.6).

Fig. 7 displays the outage status, i.e. whether we have an outage in the
selectivity results provided in Fig. 6(b) or not, versus the window number. To
create Fig. 7, we averaged the selectivity values over windows of size
$T_{s}=5$ ms. Since our simulations span $150$ ms, this results in $30$
windows. We then determined, based on the assigned threshold, whether each
window experienced an outage by checking if the selectivity value fell below
the threshold.

For $\gamma_{o}=0.6$, we observed that out of the $30$ windows, $9$ windows
experienced an outage. This results in a selectivity outage probability of
$0.3$, indicating a $30\%$ chance that the system cannot differentiate between
the rhodopsin and bacteriorhodopsin populations. Essentially, there exists a
trade-off between the selectivity threshold set for the system and the
potential for outages. A higher threshold leads to higher selectivity values
but increases the likelihood of system outages, which imposes additional power
constraints.

Figure 7: Outage versus window number when $\gamma_{o}=0.6$. A value of zero indicates an outage, while a value of 1 indicates no outage. TABLE IV: Outage state vs Window Number Threshold |  Number of Windows with Outage |  Probability of Selectivity Outage
---|---|---
0.2 |  2 |  0.067
0.4 |  3 |  0.1
0.6 |  9 |  0.3
0.8 |  20 |  0.67

###  VI-C Selectivity Optimization Results

Experimental design, careful optimization, and validation are essential to
ensure selectivity and minimize off-target effects. In fact, to design a
system with maximum selectivity, we must determine the optimal values for both
the nanoantenna transmit power and the operating frequency, as these are the
controllable parameters. Therefore, we formulate a joint optimization problem
for these parameters, as explained in Sec. V, for both fast and slow fading
scenarios.

####  VI-C1 Fast Fading

Fig. 8(a) provides a pseudo-color plot for the selectivity of the targeted
rhodopsin population in the presence of that of bacteriorhodopsin. The plot
discriminates regions that should be targeted (given in yellow) from those
that should be avoided (given in blue) in order to activate the correct
protein population. From Fig. 8(a), we conclude that by tuning the nanoantenna
to the resonant frequency of rhodopsin, i.e. $1.36$ THz, and to a Tx power of
$30~{}\mu$W, maximum selectivity is achieved.

####  VI-C2 Slow Fading

Moving to the slow fading scenario, Fig. 8(b) demonstrates a pseudo-color plot
for the selectivity of the targeted rhodopsin population. We notice here that
in comparison to Fig. 8(a), the system needs higher Tx power, around
$38~{}\mu$W, to achieve a maximum selectivity score. This power is required to
overcome the possible outages in the system that may occur in occasions of
deep fades. Through the presented results, one can determine the amount of
power that will prevent outages and concurrently prevent the activation of the
untargeted protein population. Hence, when designing biosensors for targeted
therapy applications, engineers must carefully take into account all these
tuning parameters.

Figure 8: (a) Selectivity of rhodopsin (a) under fast fading scenario. (b)
under slow fading scenario.

It should be noted that operating nanonetworks at power levels in the
microwatt range enables the creation of energy-efficient, miniaturized devices
with improved heat dissipation capabilities. This approach aligns with the
unique requirements and constraints associated with nanoscale communication
and integration, making it feasible to deploy nanonetworks in various
applications, ranging from healthcare and environmental monitoring to
industrial and scientific domains.

##  VII Conclusions

In this work, we analyze the influence of fading on THz-induced protein
interactions. We investigate how the attenuation in the transmitted THz signal
affects the ability to achieve selective protein folding behaviors in the
desired protein population. Under fast fading channel conditions, we observed
that the protein folding probability is lower than in the case of slow fading,
primarily due to the variability in the channel response

While it’s easier to track the changes imposed by the channel under slow
fading conditions, the system may still experience outages, resulting in
significant delays in receiving the nanoantenna signal and placing additional
demands on the power requirements of the system.

Through the examination of various fading scenarios achieved by varying the
velocity of RBCs, we illustrate the influence of the cellular environment on
protein dynamics. This analysis can provide insights into the effects of
molecular crowding and hindered diffusion on protein folding. By taking into
account channel variability, we can design biosensors and nanoantennas that
align with medical requirements and fully harness the potential of the intra-
body communication channel.

However, one of the current limitations of the developed model is the absence
of experimental validation. While theoretically, THz antennas could be
miniaturized due to the short wavelength of THz-EM waves, current systems
often rely on larger and more obtrusive equipment, such as tabletop
microscopes and spectroscopy systems. These factors create barriers that
constrain the potential and possibilities of this emerging field. Research
efforts are therefore needed to address these challenges and enhance the
feasibility of intra-body THz experiments for clinical applications.
Advancements in nanotechnology offer a promising solution, as they provide
diverse opportunities to enhance the performance, sensitivity, and
applicability of THz experiments across various fields.

## References

  * [1] H. Elayan, A. Eckford, and R. Adve, “The Impact of Slow Fading on THz-Induced Protein Interactions,” (accepted for publication in the IEEE International Conference on Communications, 2023).
  * [2] I. F. Akyildiz and J. M. Jornet, “The internet of nano-things,” _IEEE Wireless Communications_ , vol. 17, no. 6, pp. 58–63, 2010.
  * [3] J. M. Jornet and I. F. Akyildiz, “Graphene-based plasmonic nano-antenna for terahertz band communication in nanonetworks,” _IEEE Journal on selected areas in communications_ , vol. 31, no. 12, pp. 685–694, 2013.
  * [4] K. Yang _et al._ , “Numerical analysis and characterization of THz propagation channel for body-centric nano-communications,” _IEEE Transactions on Terahertz Science and technology_ , vol. 5, no. 3, pp. 419–426, 2015.
  * [5] H. Elayan, R. M. Shubair, J. M. Jornet, and P. Johari, “Terahertz channel model and link budget analysis for intrabody nanoscale communication,” _IEEE transactions on nanobioscience_ , vol. 16, no. 6, pp. 491–503, 2017\.
  * [6] R. Indrawijaya and T. Kürner, “Simulation of terahertz intrabody wireless nano sensor networks in the presence of noise and interference,” _12th European Conference on Antennas and Propagation_ , 2018.
  * [7] Q. H. Abbasi _et al._ , “Terahertz channel characterization inside the human skin for nano-scale body-centric networks,” _IEEE Transactions on Terahertz Science and Technology_ , vol. 6, no. 3, pp. 427–434, 2016.
  * [8] A. Sangwan _et al._ , “Increasing the communication distance between nano-biosensing implants and wearable devices,” in _2018 IEEE 19th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_. IEEE, 2018, pp. 1–5.
  * [9] J. M. Jornet and A. Sangwan, “Nanonetworking in the Terahertz Band and Beyond,” _IEEE Nanotechnology Magazine_ , 2023.
  * [10] L. Wei _et al._ , “Application of terahertz spectroscopy in biomolecule detection,” _Frontiers in Laboratory Medicine_ , vol. 2, no. 4, pp. 127–133, 2018.
  * [11] A. G. Markelz, “Terahertz dielectric sensitivity to biomolecular structure and function,” _IEEE Journal of Selected Topics in Quantum Electronics_ , vol. 14, no. 1, pp. 180–190, 2008.
  * [12] A. M. Stadler _et al._ , “Picosecond to nanosecond dynamics provide a source of conformational entropy for protein folding,” _Physical Chemistry Chemical Physics_ , vol. 18, no. 31, pp. 21 527–21 538, 2016.
  * [13] W. A. Eaton, P. A. Thompson, C.-K. Chan, S. J. Hage, and J. Hofrichter, “Fast events in protein folding,” _Structure_ , vol. 4, no. 10, pp. 1133–1139, 1996.
  * [14] K. Dave and M. Gruebele, “Fast-folding proteins under stress,” _Cellular and Molecular Life Sciences_ , vol. 72, no. 22, pp. 4273–4285, 2015\.
  * [15] H. Elayan, A. Eckford, and R. Adve, “Toward establishing molecular interfaces using terahertz radiation,” in _2022 IEEE 16th International Symposium on Medical Information and Communication Technology (ISMICT)_. IEEE, 2022, pp. 1–6.
  * [16] H. Elayan, A. W. Eckford, and R. S. Adve, “Information rates of controlled protein interactions using terahertz communication,” _IEEE Transactions on NanoBioscience_ , vol. 20, no. 1, pp. 9–19, 2020.
  * [17] ——, “Selectivity of protein interactions stimulated by Terahertz signals,” _IEEE Transactions on NanoBioscience_ , 2022.
  * [18] P. Johari and J. M. Jornet, “Nanoscale optical wireless channel model for intra-body communications: Geometrical, time, and frequency domain analyses,” _IEEE Transactions on Communications_ , vol. 66, no. 4, pp. 1579–1593, 2017.
  * [19] A. Sangwan and J. M. Jornet, “Beamforming optical antenna arrays for nano-bio sensing and actuation applications,” _Nano Communication Networks_ , vol. 29, p. 100363, 2021.
  * [20] D. K. Cheng _et al._ , _Field and wave electromagnetics_. Pearson Education India, 1989.
  * [21] Y.-S. Lee, _Principles of Terahertz Science and Technology_. Springer Science & Business Media, 2009, vol. 170\.
  * [22] D. Kim, M. Sonker, and A. Ros, “Dielectrophoresis: From molecular to micrometer-scale analytes,” _Analytical chemistry_ , vol. 91, no. 1, pp. 277–295, 2018.
  * [23] J. Knab, J.-Y. Chen, and A. Markelz, “Hydration dependence of conformational dielectric relaxation of lysozyme,” _Biophysical journal_ , vol. 90, no. 7, pp. 2576–2581, 2006.
  * [24] A. Ishimaru, _Wave propagation and scattering in random media_. Academic press New York, 1978, vol. 2.
  * [25] C. Enjamio, E. Vilar, and F. Perez-Fontan, “Rain scatter interference in mm-wave broadband fixed wireless access networks caused by a 2-D dynamic rain environment,” _IEEE transactions on wireless communications_ , vol. 6, no. 7, pp. 2497–2507, 2007.
  * [26] A. Holt, R. McGuinness, D. Charlton, P. Thompson, and M. Mehler, “The development of a model to estimate the bistatic transmission loss associated with intersystem interference,” _IEEE transactions on antennas and propagation_ , vol. 41, no. 10, pp. 1422–1431, 1993.
  * [27] H. Elayan, A. W. Eckford, and R. S. Adve, “Terahertz Intra-body Propagation through LOS and NLOS Links,” in _ICC 2022 - IEEE International Conference on Communications_ , 2022, pp. 1716–1721.
  * [28] C. E. McLaren, G. M. Brittenham, and V. Hasselblad, “Analysis of the volume of red blood cells: application of the expectation-maximization algorithm to grouped data from the doubly-truncated lognormal distribution,” _Biometrics_ , pp. 143–158, 1986.
  * [29] J. McCammon, “Protein dynamics,” _Reports on Progress in Physics_ , vol. 47, no. 1, p. 1, 1984.
  * [30] L. Cai, N. Friedman, and X. S. Xie, “Stochastic protein expression in individual cells at the single molecule level,” _Nature_ , vol. 440, no. 7082, pp. 358–362, 2006.
  * [31] V. Balakrishnan, _Elements of nonequilibrium statistical mechanics_. Ane Books, 2008, vol. 3.
  * [32] A. Carpinteri _et al._ , “Terahertz mechanical vibrations in lysozyme: Raman spectroscopy vs modal analysis,” _Journal of Molecular Structure_ , vol. 1139, pp. 222–230, 2017.
  * [33] M. M. Tirion, “Large amplitude elastic motions in proteins from a single-parameter, atomic analysis,” _Physical review letters_ , vol. 77, no. 9, p. 1905, 1996.
  * [34] F. W. Byron and R. W. Fuller, _Mathematics of classical and quantum physics_. Courier Corporation, 2012.
  * [35] D. A. Turton _et al._ , “Terahertz underdamped vibrational motion governs protein-ligand binding in solution,” _Nature communications_ , vol. 5, no. 1, pp. 1–6, 2014.
  * [36] G. Acbas _et al._ , “Optical measurements of long-range protein vibrations,” _Nature communications_ , vol. 5, no. 1, pp. 1–7, 2014.
  * [37] H. Elayan, A. W. Eckford, and R. Adve, “Enabling Protein Interactions Using Terahertz Signals for Intra-body Communication,” in _Proceedings of the 19th ACM Conference on Embedded Networked Sensor Systems_ , 2021, pp. 603–609.
  * [38] A. V. Finkelstein _et al._ , “Why do protein architectures have Boltzmann-like statistics?” _Proteins: Structure, Function, and Bioinformatics_ , vol. 23, no. 2, pp. 142–150, 1995.
  * [39] F. Sachs and H. Lecar, “Stochastic models for mechanical transduction,” _Biophysical journal_ , vol. 59, no. 5, pp. 1143–1145, 1991.
  * [40] C. B. Anfinsen, “Principles that govern the folding of protein chains,” _Science_ , vol. 181, no. 4096, pp. 223–230, 1973.
  * [41] D. Tse and P. Viswanath, _Fundamentals of wireless communication_. Cambridge university press, 2005.
  * [42] A. Cegielski, _Iterative methods for fixed point problems in Hilbert spaces_. Springer, 2012, vol. 2057.
  * [43] R. Balu _et al._ , “Terahertz spectroscopy of bacteriorhodopsin and rhodopsin: similarities and differences,” _Biophysical journal_ , vol. 94, no. 8, pp. 3217–3226, 2008.
  * [44] K. T. Sapra, P. S.-H. Park, K. Palczewski, and D. J. Muller, “Mechanical Properties of Bovine Rhodopsin and Bacteriorhodopsin: Possible Roles in Folding and Function,” _Langmuir_ , vol. 24, no. 4, pp. 1330–1337, 2008\.
  * [45] C. J. Benham _et al._ , _Mathematics of DNA structure, function and interactions_. Springer Science & Business Media, 2009, vol. 150.
  * [46] C. B. Reid _et al._ , “Terahertz time-domain spectroscopy of human blood,” _IEEE journal of biomedical and health informatics_ , vol. 17, no. 4, pp. 774–778, 2013.
  * [47] B. Manavalan, K. Kuwajima, and J. Lee, “PFDB: a standardized protein folding database with temperature correction,” _Scientific reports_ , vol. 9, no. 1, pp. 1–9, 2019.
  * [48] E. Chaigneau _et al._ , “Unbiased analysis method for measurement of red blood cell size and velocity with laser scanning microscopy,” _Frontiers in neuroscience_ , p. 644, 2019.


except in the extended scenarios assuming a cosmological constant, where for
CMB+Pantheon case $S_{8}$ values are slightly close to it. For all the other
cases, for these two datasets the tension in the $S_{8}$ parameter increases
significantly. This is because of the indication for a closed Universe coming
from these two datasets when it is not assumed a cosmological constant. On the
other hand, when we look at the estimations of $S_{8}$ by CMB+BAO and
CMB+BAO+Pantheon, we see that these estimations are almost similar to the
Planck values assuming a vanilla $\Lambda$CDM model. This indicates that
within such extended scenarios, the $S_{8}$ tension increases when there is an
indication for a closed Universe.

####  V.1.4 $M_{\nu}$

From all the analyses, we find no evidence of a total neutrino mass. In
particular, only in a few cases, namely, $\Lambda$CDM + $\Omega_{K}$ \+
$M_{\nu}$ \+ $N_{\rm eff}$, $w$CDM + $\Omega_{K}$ \+ $M_{\nu}$,
$w_{0}w_{a}$CDM + $\Omega_{K}$ \+ $M_{\nu}$, PEDE + $\Omega_{K}$ \+ $M_{\nu}$
\+ $N_{\rm eff}$, for CMB alone, we find an indication of a total neutrino
mass at 68% CL but this evidence disappears when other external and
complementary probes are added to CMB. However, in order to realize the bounds
on $M_{\nu}$ in all these scenarios, in Figure 20, we have presented a
graphical presentation of the the upper limits of $M_{\nu}$ at 95% CL for our
extended cosmological models using CMB, CMB+BAO, CMB+Pantheon and
CMB+BAO+Pantheon.

Figure 20: Whisker plot with 95% CL upper bounds on $M_{\nu}$ [eV] obtained in
various cosmological scenarios considered in this article has been displayed
for several observational datasets, namely, CMB, CMB+BAO, CMB+Pantheon and
CMB+BAO+Pantheon.

####  V.1.5 $N_{\rm eff}$

We provide a summary picture of $N_{\rm eff}$ in Figure 21 obtained in all our
extended cosmological scenarios considered in this work for CMB, CMB+BAO,
CMB+Pantheon, and CMB+BAO+Pantheon. We find that irrespective of all the
cosmological models considered in this work, the estimated values of $N_{\rm
eff}$ for CMB and CMB+Pantheon are almost same with that of the standard value
$N_{\rm eff}=3.044$ and thus these estimations are consistent with this
standard value. For the remaining two cases, namely, CMB+BAO and
CMB+BAO+Pantheon, the values of $N_{\rm eff}$ are slightly lower ($N_{\rm
eff}<3$) but within 68% CL, they are also consistent with the standard value.

Figure 21: Whisker plot with 68% CL constraints on $N_{\rm eff}$ obtained in
various cosmological scenarios considered in this article has been displayed
for several observational datasets, namely, CMB, CMB+BAO, CMB+Pantheon and
CMB+BAO+Pantheon. The magenta vertical line corresponds to the standard value
$N_{\rm eff}=3.044$.

####  V.1.6 $w_{0}$ and $w_{a}$

Here we summarize the main results on the dark energy equation of state
parameter $w$ of the $w$CDM model and also on the $(w_{0},w_{a})$ parameters
of the $w_{0}w_{a}$CDM model. In Figure 22, we show the whisker plots with 68%
CL, and 95% upper limits, on $w$ and $w_{a}$ for the $w$CDM and
$w_{0}w_{a}$CDM models for all the observational datasets.

One can clearly see that in the $w$CDM model, only the CMB+Pantheon dataset
gives a strong indication of a phantom DE (i.e. $w<-1$) at more than 95% CL,
only when the equation of state of the DE is constant and not dynamical. For
the remaining cases, $w=-1$ is allowed within 68% CL.

However, for the $w_{0}w_{a}$CDM model in which the dark energy equation of
state is dynamical, CMB+BAO gives an indication for a quintessence DE (i.e.
$w_{0}>-1$), while all the other dataset combinations are in agreement with a
cosmological constant777Note here that $w_{0}$ remains unconstrained in the
scenarios $w_{0}w_{a}$CDM + $\Omega_{K}$ \+ $M_{\nu}$, $w_{0}w_{a}$CDM +
$\Omega_{K}$ \+ $M_{\nu}$ \+ $N_{\rm eff}$ and for the scenario
$w_{0}w_{a}$CDM + $\Omega_{K}$ \+ $N_{\rm eff}$, even though the 68% CL upper
bound on $w_{0}$ is available but it becomes unconstrained at 95% CL. .
Regarding the dynamical nature of the DE quantified through the free parameter
$w_{a}$, we find that CMB alone and CMB+Pantheon are very consistent with
$w_{a}=0$ (except in the model $w_{0}w_{a}$CDM + $\Omega_{K}$ in which
$w_{a}\neq 0$ at more than 68% CL), however, when BAO data are included to
CMB, an indication for a $w_{a}<0$ at more than $2\sigma$ is found. While this
evidence is mildly diluted (but remains true at slightly more than 68% CL)
when the full combination CMB+BAO+Pantheon is explored. In other words, for
this $w_{0}w_{a}$CDM model, CMB+BAO is preferring a quintessence dynamical DE
at more than 95% CL.

Figure 22: Whisker plot with 68% CL constraints and 95% CL limits on $w_{0}$
and $w_{a}$ obtained in various cosmological scenarios considered in this
article has been displayed for several observational datasets, namely, CMB,
CMB+BAO, CMB+Pantheon and CMB+BAO+Pantheon. The magenta vertical lines in the
two panels represent to the standard cosmological constant value $w_{0}=-1$
and to the non dynamical $w_{a}=0$ limit, respectively. The missing blue lines
(Planck only results) in the left panel are indicating that $w$ is
unconstrained for those models.

###  V.2 Concluding remarks and the road ahead

There is no doubt that without observational data, modern cosmology is
incomplete. From the detection of the cosmic microwave background radiation to
the discovery of the late-time cosmic acceleration, we have witnessed how the
observational data have been crucially important in the understanding of our
Universe. Whenever we talk about cosmology, the six parameter $\Lambda$CDM
cosmological model under the assumptions of a positive cosmological constant,
a pressure-less DM and a flat Universe, naturally enters into the picture
because of its super agreement with a series of astronomical data. However,
the true nature of DE and DM within this simplest cosmological scenario has
remained one of the longstanding debates in cosmology. In addition to that,
the assumption of a flat Universe has created a further debate in the
scientific community. A series of recent articles have argued that the present
observational data hint towards an evidence of a closed Universe Aghanim  _et
al._ (2020a); Handley (2019); Di Valentino  _et al._ (2019, 2020b); Glanville
_et al._ (2022). This evidence has sparked the entire scientific community
because any evidence of a closed Universe increases the tensions in the key
cosmological parameters and hence the flat $\Lambda$CDM cosmology is further
challenged. Generally, the curvature of the Universe is “assumed” to be zero
when the cosmological models are analyzed. This is motivated because: (i) the
inflationary paradigm has been tremendously successful in demanding a
spatially flat Universe, (ii) the observational data in the past preferred a
small curvature of the Universe which is very close to zero, and further (iii)
the inclusion of extra parameters, in terms of the curvature parameter for
instance, may increase the degeneracies in the parameters. However, over the
last several years, the sensitivity of the astronomical data has improved a
lot, and given such improvements, it is very natural to allow the
observational data to decide the nature of the curvature parameter. That
means, the consideration of a flat Universe (i.e. $\Omega_{K}=0$) while
analyzing the cosmological models is not a very realistic assumption. We
further comment that forcing $\Omega_{K}=0$ could bias the results and the
intrinsic nature of the cosmological models may remain hidden forever Anselmi
_et al._ (2022), and at this stage where we have a large amount of
astronomical data from various sources, having $\Omega_{K}=0$ is not expected.

Thus, in the present article we investigated the effects of curvature on the
cosmological models and their associated (free and derived) parameters. We
considered some classical and well known cosmological models and their several
extensions, including variations in the DE and neutrino sectors, in presence
of the curvature of our Universe. Our analyses clearly indicate that it is
better not to fix $\Omega_{K}=0$, to analyze the cosmological models. We trust
that the upcoming observational data from various astronomical surveys (see
for instance, Refs. Abell  _et al._ (2009); Laureijs  _et al._ (2011); Levi
_et al._ (2013); Ade  _et al._ (2019)) will be highly promising in this
direction which will reveal the intrinsic nature of the dark sector as well as
the curvature of the Universe.

###### Acknowledgements.

The authors thank the referee for some useful comments which improved the
quality of discussion of the article. WY was supported by the National Natural
Science Foundation of China under Grants No. 12175096 and No. 11705079, and
Liaoning Revitalization Talents Program under Grant no. XLYC1907098. WG and AM
are supported by “Theoretical Astroparticle Physics” (TAsP), iniziativa
specifica INFN. SP acknowledges the financial support from the Department of
Science and Technology (DST), Govt. of India under the Scheme “Fund for
Improvement of S&T Infrastructure (FIST)” (File No. SR/FST/MS-I/2019/41). EDV
is supported by a Royal Society Dorothy Hodgkin Research Fellowship.

## References

  * Di Valentino  _et al._ (2020a) E. Di Valentino _et al._ , (2020a), arXiv:2008.11286 [astro-ph.CO] .
  * Gaztanaga  _et al._ (2009) E. Gaztanaga, R. Miquel, and E. Sanchez, Phys. Rev. Lett. 103, 091302 (2009), arXiv:0808.1921 [astro-ph] .
  * Mortonson (2009) M. J. Mortonson, Phys. Rev. D 80, 123504 (2009), arXiv:0908.0346 [astro-ph.CO] .
  * Suyu  _et al._ (2014) S. H. Suyu _et al._ , Astrophys. J. Lett. 788, L35 (2014), arXiv:1306.4732 [astro-ph.CO] .
  * L’Huillier and Shafieloo (2017) B. L’Huillier and A. Shafieloo, JCAP 01, 015 (2017), arXiv:1606.06832 [astro-ph.CO] .
  * Chudaykin  _et al._ (2021) A. Chudaykin, K. Dolgikh, and M. M. Ivanov, Phys. Rev. D 103, 023507 (2021), arXiv:2009.10106 [astro-ph.CO] .
  * Acquaviva  _et al._ (2021) G. Acquaviva, O. Akarsu, N. Katirci, and J. A. Vazquez, Phys. Rev. D 104, 023505 (2021), arXiv:2104.02623 [astro-ph.CO] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020a), arXiv:1807.06209 [astro-ph.CO] .
  * Handley (2019) W. Handley, arXiv:1908.09139 (2019), arXiv:1908.09139 [astro-ph.CO] .
  * Di Valentino  _et al._ (2019) E. Di Valentino, A. Melchiorri, and J. Silk, Nature Astron. 4, 196 (2019), arXiv:1911.02087 [astro-ph.CO] .
  * Di Valentino  _et al._ (2020b) E. Di Valentino, A. Melchiorri, and J. Silk, arXiv:2003.04935 (2020b), arXiv:2003.04935 [astro-ph.CO] .
  * Semenaite  _et al._ (2022) A. Semenaite, A. G. Sánchez, A. Pezzotta, J. Hou, A. Eggemeier, M. Crocce, C. Zhao, J. R. Brownstein, G. Rossi, and D. P. Schneider, (2022), arXiv:2210.07304 [astro-ph.CO] .
  * Calabrese  _et al._ (2008) E. Calabrese, A. Slosar, A. Melchiorri, G. F. Smoot, and O. Zahn, Phys. Rev. D77, 123531 (2008), arXiv:0803.2309 [astro-ph] .
  * Di Valentino  _et al._ (2020c) E. Di Valentino, A. Melchiorri, and J. Silk, JCAP 01, 013 (2020c), arXiv:1908.01391 [astro-ph.CO] .
  * Efstathiou and Gratton (2020) G. Efstathiou and S. Gratton,  (2020), 10.1093/mnrasl/slaa093, arXiv:2002.06892 [astro-ph.CO] .
  * Efstathiou and Gratton (2019) G. Efstathiou and S. Gratton, (2019), arXiv:1910.00483 [astro-ph.CO] .
  * Rosenberg  _et al._ (2022) E. Rosenberg, S. Gratton, and G. Efstathiou, (2022), arXiv:2205.10869 [astro-ph.CO] .
  * Glanville  _et al._ (2022) A. Glanville, C. Howlett, and T. M. Davis, (2022), arXiv:2205.05892 [astro-ph.CO] .
  * Anselmi  _et al._ (2022) S. Anselmi, M. F. Carney, J. T. Giblin, S. Kumar, J. B. Mertens, M. ODwyer, G. D. Starkman, and C. Tian, (2022), arXiv:2207.06547 [astro-ph.CO] .
  * Dossett and Ishak (2012) J. Dossett and M. Ishak, Phys. Rev. D 86, 103008 (2012), arXiv:1205.2422 [astro-ph.CO] .
  * Abdalla  _et al._ (2022) E. Abdalla _et al._ , JHEAp 34, 49 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Cao and Ratra (2022) S. Cao and B. Ratra, Mon. Not. Roy. Astron. Soc. 513, 5686 (2022), arXiv:2203.10825 [astro-ph.CO] .
  * Qi  _et al._ (2022) J.-Z. Qi, Y. Cui, W.-H. Hu, J.-F. Zhang, J.-L. Cui, and X. Zhang, Phys. Rev. D 106, 023520 (2022), arXiv:2202.01396 [astro-ph.CO] .
  * Cao  _et al._ (2022a) S. Cao, N. Khadka, and B. Ratra, Mon. Not. Roy. Astron. Soc. 510, 2928 (2022a), arXiv:2110.14840 [astro-ph.CO] .
  * Cao  _et al._ (2022b) S. Cao, J. Ryan, and B. Ratra, Mon. Not. Roy. Astron. Soc. 509, 4745 (2022b), arXiv:2109.01987 [astro-ph.CO] .
  * Cao  _et al._ (2022c) M.-D. Cao, J. Zheng, J.-Z. Qi, X. Zhang, and Z.-H. Zhu, Astrophys. J. 934, 108 (2022c), arXiv:2112.14564 [astro-ph.CO] .
  * Cao  _et al._ (2021) S. Cao, J. Ryan, and B. Ratra, Mon. Not. Roy. Astron. Soc. 504, 300 (2021), arXiv:2101.08817 [astro-ph.CO] .
  * Qi  _et al._ (2021) J.-Z. Qi, J.-W. Zhao, S. Cao, M. Biesiada, and Y. Liu, Mon. Not. Roy. Astron. Soc. 503, 2179 (2021), arXiv:2011.00713 [astro-ph.CO] .
  * Cao  _et al._ (2020) S. Cao, J. Ryan, and B. Ratra, Mon. Not. Roy. Astron. Soc. 497, 3191 (2020), arXiv:2005.12617 [astro-ph.CO] .
  * Wang  _et al._ (2020) B. Wang, J.-Z. Qi, J.-F. Zhang, and X. Zhang, Astrophys. J. 898, 100 (2020), arXiv:1910.12173 [astro-ph.CO] .
  * Riess  _et al._ (2022a) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022a), arXiv:2112.04510 [astro-ph.CO] .
  * Riess  _et al._ (2022b) A. G. Riess, L. Breuval, W. Yuan, S. Casertano, L. M. ~Macri, D. Scolnic, T. Cantat-Gaudin, R. I. Anderson, and M. C. Reyes, (2022b), arXiv:2208.01045 [astro-ph.CO] .
  * Di Valentino  _et al._ (2020d) E. Di Valentino _et al._ , arXiv:2008.11284 (2020d), arXiv:2008.11284 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021a) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021a), arXiv:2103.01183 [astro-ph.CO] .
  * Di Valentino  _et al._ (2017a) E. Di Valentino, A. Melchiorri, and O. Mena, Phys. Rev. D 96, 043503 (2017a), arXiv:1704.08342 [astro-ph.CO] .
  * Kumar and Nunes (2017) S. Kumar and R. C. Nunes, Phys. Rev. D96, 103511 (2017), arXiv:1702.02143 [astro-ph.CO] .
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A. G. Riess, in  _Nature Astronomy 2019_ (2019) arXiv:1907.10625 [astro-ph.CO] .
  * Knox and Millea (2020) L. Knox and M. Millea, Phys. Rev. D 101, 043533 (2020), arXiv:1908.03663 [astro-ph.CO] .
  * Jedamzik  _et al._ (2021) K. Jedamzik, L. Pogosian, and G.-B. Zhao, Commun. in Phys. 4, 123 (2021), arXiv:2010.04158 [astro-ph.CO] .
  * Di Valentino  _et al._ (2018) E. Di Valentino, C. Bøehm, E. Hivon, and F. c. R. Bouchet, Phys. Rev. D 97, 043513 (2018), arXiv:1710.02559 [astro-ph.CO] .
  * Yang  _et al._ (2018a) W. Yang, S. Pan, E. Di Valentino, R. C. Nunes, S. Vagnozzi, and D. F. Mota, JCAP 09, 019 (2018a), arXiv:1805.08252 [astro-ph.CO] .
  * Yang  _et al._ (2018b) W. Yang, A. Mukherjee, E. Di Valentino, and S. Pan, Phys. Rev. D98, 123527 (2018b), arXiv:1809.06883 [astro-ph.CO] .
  * Pan  _et al._ (2019a) S. Pan, W. Yang, C. Singha, and E. N. Saridakis, Phys. Rev. D100, 083539 (2019a), arXiv:1903.10969 [astro-ph.CO] .
  * Yang  _et al._ (2020a) W. Yang, S. Pan, R. C. Nunes, and D. F. Mota, JCAP 04, 008 (2020a), arXiv:1910.08821 [astro-ph.CO] .
  * Pan  _et al._ (2019b) S. Pan, W. Yang, E. Di Valentino, E. N. Saridakis, and S. Chakraborty, Phys. Rev. D 100, 103520 (2019b), arXiv:1907.07540 [astro-ph.CO] .
  * Poulin  _et al._ (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Yang  _et al._ (2019a) W. Yang, S. Pan, E. Di Valentino, E. N. Saridakis, and S. Chakraborty, Phys. Rev. D99, 043543 (2019a), arXiv:1810.05141 [astro-ph.CO] .
  * Pan  _et al._ (2020a) S. Pan, W. Yang, and A. Paliathanasis, Mon. Not. Roy. Astron. Soc. 493, 3114 (2020a), arXiv:2002.03408 [astro-ph.CO] .
  * Di Valentino  _et al._ (2020e) E. Di Valentino, A. Melchiorri, O. Mena, and S. Vagnozzi, Phys. Dark Univ. 30, 100666 (2020e), arXiv:1908.04281 [astro-ph.CO] .
  * Di Valentino  _et al._ (2020f) E. Di Valentino, A. Melchiorri, O. Mena, and S. Vagnozzi, Phys. Rev. D 101, 063502 (2020f), arXiv:1910.09853 [astro-ph.CO] .
  * Yao and Meng (2020) Y. Yao and X. Meng, (2020), arXiv:2011.09160 [astro-ph.CO] .
  * Lucca and Hooper (2020) M. Lucca and D. C. Hooper, Phys. Rev. D 102, 123502 (2020), arXiv:2002.06127 [astro-ph.CO] .
  * Blinov  _et al._ (2020) N. Blinov, C. Keith, and D. Hooper, JCAP 06, 005 (2020), arXiv:2004.06114 [astro-ph.CO] .
  * Anchordoqui  _et al._ (2021) L. A. Anchordoqui, E. Di Valentino, S. Pan, and W. Yang, JHEAp 32, 28 (2021), arXiv:2107.13932 [astro-ph.CO] .
  * Karwal  _et al._ (2022) T. Karwal, M. Raveri, B. Jain, J. Khoury, and M. Trodden, Phys. Rev. D 105, 063535 (2022), arXiv:2106.13290 [astro-ph.CO] .
  * Freese and Winkler (2021) K. Freese and M. W. Winkler, Phys. Rev. D 104, 083533 (2021), arXiv:2102.13655 [astro-ph.CO] .
  * Perivolaropoulos and Skara (2022) L. Perivolaropoulos and F. Skara, New Astron. Rev. 95, 101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Schöneberg  _et al._ (2022) N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, Phys. Rept. 984, 1 (2022), arXiv:2107.10291 [astro-ph.CO] .
  * Reeves  _et al._ (2022) A. Reeves, L. Herold, S. Vagnozzi, B. D. Sherwin, and E. G. M. Ferreira, (2022), arXiv:2207.01501 [astro-ph.CO] .
  * Colgáin  _et al._ (2022a) E. O. Colgáin, M. M. Sheikh-Jabbari, R. Solomon, G. Bargiacchi, S. Capozziello, M. G. Dainotti, and D. Stojkovic, Phys. Rev. D 106, L041301 (2022a), arXiv:2203.10558 [astro-ph.CO] .
  * Naidoo  _et al._ (2022) K. Naidoo, M. Jaber, W. A. Hellwing, and M. Bilicki, (2022), arXiv:2209.08102 [astro-ph.CO] .
  * Cruz  _et al._ (2022) J. S. Cruz, F. Niedermann, and M. S. Sloth, (2022), arXiv:2209.02708 [astro-ph.CO] .
  * Escudero  _et al._ (2022) H. G. Escudero, J.-L. Kuo, R. E. Keeley, and K. N. Abazajian, (2022), arXiv:2208.14435 [astro-ph.CO] .
  * Gómez-Valent  _et al._ (2022) A. Gómez-Valent, Z. Zheng, L. Amendola, C. Wetterich, and V. Pettorino, (2022), arXiv:2207.14487 [astro-ph.CO] .
  * Colgáin  _et al._ (2022b) E. O. Colgáin, M. M. Sheikh-Jabbari, R. Solomon, M. G. Dainotti, and D. Stojkovic, (2022b), arXiv:2206.11447 [astro-ph.CO] .
  * Yang  _et al._ (2019b) W. Yang, S. Pan, A. Paliathanasis, S. Ghosh, and Y. Wu, Mon. Not. Roy. Astron. Soc. 490, 2071 (2019b), arXiv:1904.10436 [gr-qc] .
  * Yang  _et al._ (2020b) W. Yang, E. Di Valentino, S. Pan, S. Basilakos, and A. Paliathanasis, Phys. Rev. D 102, 063503 (2020b), arXiv:2001.04307 [astro-ph.CO] .
  * Di Valentino  _et al._ (2022a) E. Di Valentino, S. Gariazzo, C. Giunti, O. Mena, S. Pan, and W. Yang, Phys. Rev. D 105, 103511 (2022a), arXiv:2110.03990 [astro-ph.CO] .
  * Yang  _et al._ (2022) W. Yang, S. Pan, O. Mena, and E. Di Valentino, (2022), arXiv:2209.14816 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021b) E. Di Valentino _et al._ , Astropart. Phys. 131, 102604 (2021b), arXiv:2008.11285 [astro-ph.CO] .
  * Heymans  _et al._ (2021) C. Heymans _et al._ , Astron. Astrophys. 646, A140 (2021), arXiv:2007.15632 [astro-ph.CO] .
  * Tröster  _et al._ (2021) T. Tröster _et al._ (KiDS), Astron. Astrophys. 649, A88 (2021), arXiv:2010.16416 [astro-ph.CO] .
  * Secco  _et al._ (2022) L. F. Secco _et al._ (DES), Phys. Rev. D 105, 023515 (2022), arXiv:2105.13544 [astro-ph.CO] .
  * Abbott  _et al._ (2022) T. M. C. Abbott _et al._ (DES), (2022), arXiv:2207.05766 [astro-ph.CO] .
  * Franco Abellán  _et al._ (2022) G. Franco Abellán, R. Murgia, V. Poulin, and J. Lavalle, Phys. Rev. D 105, 063525 (2022), arXiv:2008.09615 [astro-ph.CO] .
  * Franco Abellán  _et al._ (2021) G. Franco Abellán, R. Murgia, and V. Poulin, Phys. Rev. D 104, 123533 (2021), arXiv:2102.12498 [astro-ph.CO] .
  * Lucca (2021) M. Lucca, Phys. Dark Univ. 34, 100899 (2021), arXiv:2105.09249 [astro-ph.CO] .
  * de Araujo  _et al._ (2021) J. C. N. de Araujo, A. De Felice, S. Kumar, and R. C. Nunes, Phys. Rev. D 104, 104057 (2021), arXiv:2106.09595 [astro-ph.CO] .
  * Clark  _et al._ (2021) S. J. Clark, K. Vattis, J. Fan, and S. M. Koushiappas, (2021), arXiv:2110.09562 [astro-ph.CO] .
  * Beltrán Jiménez  _et al._ (2021) J. Beltrán Jiménez, D. Bettoni, D. Figueruelo, F. A. Teppa Pannia, and S. Tsujikawa, Phys. Rev. D 104, 103503 (2021), arXiv:2106.11222 [astro-ph.CO] .
  * Heimersheim  _et al._ (2020) S. Heimersheim, N. Schöneberg, D. C. Hooper, and J. Lesgourgues, JCAP 12, 016 (2020), arXiv:2008.08486 [astro-ph.CO] .
  * Poulin  _et al._ (2022) V. Poulin, J. L. Bernal, E. Kovetz, and M. Kamionkowski, (2022), arXiv:2209.06217 [astro-ph.CO] .
  * Amon and Efstathiou (2022) A. Amon and G. Efstathiou,  (2022), 10.1093/mnras/stac2429, arXiv:2206.11794 [astro-ph.CO] .
  * Chevallier and Polarski (2001) M. Chevallier and D. Polarski, Int. J. Mod. Phys. D10, 213 (2001), arXiv:gr-qc/0009008 [gr-qc] .
  * Linder (2003) E. V. Linder, Phys. Rev. Lett. 90, 091301 (2003), arXiv:astro-ph/0208512 [astro-ph] .
  * Li and Shafieloo (2019) X. Li and A. Shafieloo, Astrophys. J. Lett. 883, L3 (2019), arXiv:1906.08275 [astro-ph.CO] .
  * Pan  _et al._ (2020b) S. Pan, W. Yang, E. Di Valentino, A. Shafieloo, and S. Chakraborty, JCAP 06, 062 (2020b), arXiv:1907.12551 [astro-ph.CO] .
  * Perenon  _et al._ (2022) L. Perenon, M. Martinelli, R. Maartens, S. Camera, and C. Clarkson, Phys. Dark Univ. 37, 101119 (2022), arXiv:2206.12375 [astro-ph.CO] .
  * Yang  _et al._ (2021a) W. Yang, E. Di Valentino, S. Pan, A. Shafieloo, and X. Li, Phys. Rev. D 104, 063521 (2021a), arXiv:2103.03815 [astro-ph.CO] .
  * Yang  _et al._ (2021b) W. Yang, E. Di Valentino, S. Pan, Y. Wu, and J. Lu, Mon. Not. Roy. Astron. Soc. 501, 5845 (2021b), arXiv:2101.02168 [astro-ph.CO] .
  * Menci  _et al._ (2020) N. Menci _et al._ , Astrophys. J. 900, 108 (2020), arXiv:2007.12453 [astro-ph.CO] .
  * Zhao  _et al._ (2020) Z.-W. Zhao, Z.-X. Li, J.-Z. Qi, H. Gao, J.-F. Zhang, and X. Zhang, Astrophys. J. 903, 83 (2020), arXiv:2006.01450 [astro-ph.CO] .
  * Zimdahl  _et al._ (2020) W. Zimdahl, J. Fabris, H. Velten, and R. Herrera, Phys. Dark Univ. 30, 100681 (2020), arXiv:1911.12084 [astro-ph.CO] .
  * Du  _et al._ (2019) M. Du, W. Yang, L. Xu, S. Pan, and D. F. Mota, Phys. Rev. D 100, 043535 (2019), arXiv:1812.01440 [astro-ph.CO] .
  * Vagnozzi  _et al._ (2018) S. Vagnozzi, S. Dhawan, M. Gerbino, K. Freese, A. Goobar, and O. Mena, Phys. Rev. D 98, 083501 (2018), arXiv:1801.08553 [astro-ph.CO] .
  * Li  _et al._ (2018) X.-D. Li, C. G. Sabiu, C. Park, Y. Wang, G.-b. Zhao, H. Park, A. Shafieloo, J. Kim, and S. E. Hong, Astrophys. J. 856, 88 (2018), arXiv:1803.01851 [astro-ph.CO] .
  * Di Valentino  _et al._ (2017b) E. Di Valentino, A. Melchiorri, E. V. Linder, and J. Silk, Phys. Rev. D 96, 023523 (2017b), arXiv:1704.00762 [astro-ph.CO] .
  * Pan  _et al._ (2018) S. Pan, E. N. Saridakis, and W. Yang, Phys. Rev. D 98, 063510 (2018), arXiv:1712.05746 [astro-ph.CO] .
  * Yang  _et al._ (2018c) W. Yang, S. Pan, and A. Paliathanasis, Mon. Not. Roy. Astron. Soc. 475, 2605 (2018c), arXiv:1708.01717 [gr-qc] .
  * Yang  _et al._ (2017) W. Yang, R. C. Nunes, S. Pan, and D. F. Mota, Phys. Rev. D 95, 103522 (2017), arXiv:1703.02556 [astro-ph.CO] .
  * Zhao  _et al._ (2017) G.-B. Zhao _et al._ , Nature Astron. 1, 627 (2017), arXiv:1701.08165 [astro-ph.CO] .
  * Ma and Zhang (2011) J.-Z. Ma and X. Zhang, Phys. Lett. B 699, 233 (2011), arXiv:1102.2671 [astro-ph.CO] .
  * Di Valentino  _et al._ (2022b) E. Di Valentino, W. Giarè, A. Melchiorri, and J. Silk, (2022b), arXiv:2209.12872 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A5 (2020b), arXiv:1907.12875 [astro-ph.CO] .
  * Beutler  _et al._ (2011) F. Beutler, C. Blake, M. Colless, D. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. Roy. Astron. Soc. 416, 3017 (2011), arXiv:1106.3366 [astro-ph.CO] .
  * Ross  _et al._ (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Scolnic  _et al._ (2018) D. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Lewis and Bridle (2002) A. Lewis and S. Bridle, Phys. Rev. D 66, 103511 (2002), arXiv:astro-ph/0205436 .
  * Lewis  _et al._ (2000) A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), arXiv:astro-ph/9911177 .
  * Gelman and Rubin (1992) A. Gelman and D. B. Rubin, Statist. Sci. 7, 457 (1992).
  * Vagnozzi  _et al._ (2021) S. Vagnozzi, E. Di Valentino, S. Gariazzo, A. Melchiorri, O. Mena, and J. Silk, Phys. Dark Univ. 33, 100851 (2021), arXiv:2010.02230 [astro-ph.CO] .
  * Mangano  _et al._ (2005) G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico, Nucl. Phys. B 729, 221 (2005), arXiv:hep-ph/0506164 .
  * de Salas and Pastor (2016) P. F. de Salas and S. Pastor, JCAP 07, 051 (2016), arXiv:1606.06986 [hep-ph] .
  * Akita and Yamaguchi (2020) K. Akita and M. Yamaguchi, JCAP 08, 012 (2020), arXiv:2005.07047 [hep-ph] .
  * Froustey  _et al._ (2020) J. Froustey, C. Pitrou, and M. C. Volpe, JCAP 12, 015 (2020), arXiv:2008.01074 [hep-ph] .
  * Bennett  _et al._ (2020) J. J. Bennett, G. Buldgen, P. F. de Salas, M. Drewes, S. Gariazzo, S. Pastor, and Y. Y. Y. Wong, “Towards a precision calculation of $N_{\rm eff}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED,” (2020).
  * Archidiacono  _et al._ (2011) M. Archidiacono, E. Calabrese, and A. Melchiorri, Phys. Rev. D 84, 123008 (2011), arXiv:1109.2767 [astro-ph.CO] .
  * Di Valentino  _et al._ (2012) E. Di Valentino, M. Lattanzi, G. Mangano, A. Melchiorri, and P. Serpico, Phys. Rev. D 85, 043511 (2012), arXiv:1111.3810 [astro-ph.CO] .
  * Di Valentino  _et al._ (2013) E. Di Valentino, A. Melchiorri, and O. Mena, JCAP 11, 018 (2013), arXiv:1304.5981 [astro-ph.CO] .
  * Di Valentino  _et al._ (2016) E. Di Valentino, E. Giusarma, M. Lattanzi, O. Mena, A. Melchiorri, and J. Silk, Phys. Lett. B 752, 182 (2016), arXiv:1507.08665 [astro-ph.CO] .
  * Giarè  _et al._ (2021a) W. Giarè, E. Di Valentino, A. Melchiorri, and O. Mena, Mon. Not. Roy. Astron. Soc. 505, 2703 (2021a), arXiv:2011.14704 [astro-ph.CO] .
  * Giarè  _et al._ (2021b) W. Giarè, F. Renzi, A. Melchiorri, O. Mena, and E. Di Valentino, Mon. Not. Roy. Astron. Soc. 511, 1373 (2021b), arXiv:2110.00340 [astro-ph.CO] .
  * D’Eramo  _et al._ (2022) F. D’Eramo, E. Di Valentino, W. Giarè, F. Hajkarim, A. Melchiorri, O. Mena, F. Renzi, and S. Yun, JCAP 09, 022 (2022), arXiv:2205.07849 [astro-ph.CO] .
  * Baumann  _et al._ (2016) D. Baumann, D. Green, and B. Wallisch, Phys. Rev. Lett. 117, 171301 (2016), arXiv:1604.08614 [astro-ph.CO] .
  * Gariazzo  _et al._ (2016) S. Gariazzo, C. Giunti, M. Laveder, Y. F. Li, and E. M. Zavanin, J. Phys. G 43, 033001 (2016), arXiv:1507.08204 [hep-ph] .
  * Archidiacono and Gariazzo (2022) M. Archidiacono and S. Gariazzo, “Two sides of the same coin: sterile neutrinos and dark radiation. Status and perspectives,” (2022), arXiv:2201.10319 [hep-ph] .
  * An  _et al._ (2022) R. An, V. Gluscevic, E. Calabrese, and J. C. Hill, “What does cosmology tell us about the mass of thermal-relic dark matter?” (2022), arXiv:2202.03515 [astro-ph.CO] .
  * Abitbol  _et al._ (2019) M. H. Abitbol _et al._ (Simons Observatory), Bull. Am. Astron. Soc. 51, 147 (2019), arXiv:1907.08284 [astro-ph.IM] .
  * Di Valentino  _et al._ (2022c) E. Di Valentino, W. Giarè, A. Melchiorri, and J. Silk, (2022c), arXiv:2209.14054 [astro-ph.CO] .
  * Keeley and Shafieloo (2022) R. E. Keeley and A. Shafieloo, (2022), arXiv:2206.08440 [astro-ph.CO] .
  * Arendse  _et al._ (2020) N. Arendse _et al._ , Astron. Astrophys. 639, A57 (2020), arXiv:1909.07986 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021c) E. Di Valentino, A. Mukherjee, and A. A. Sen, Entropy 23, 404 (2021c), arXiv:2005.12587 [astro-ph.CO] .
  * Chudaykin  _et al._ (2022) A. Chudaykin, D. Gorbunov, and N. Nedelko, (2022), arXiv:2203.03666 [astro-ph.CO] .
  * Abell  _et al._ (2009) P. A. Abell _et al._ (LSST Science, LSST Project), (2009), arXiv:0912.0201 [astro-ph.IM] .
  * Laureijs  _et al._ (2011) R. Laureijs _et al._ (EUCLID), (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Levi  _et al._ (2013) M. Levi _et al._ (DESI), (2013), arXiv:1308.0847 [astro-ph.CO] .
  * Ade  _et al._ (2019) P. Ade _et al._ (Simons Observatory), JCAP 02, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .