Terahertz Induced Protein Interactions in a Random Medium

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]

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

Here, $\gamma$ is given by

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

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]

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].

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

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

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]

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

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

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]

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}$ 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

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

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

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

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

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

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.

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]

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

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

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

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

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

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

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

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

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.

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

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

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.

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.

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].

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).

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.

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.

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.

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.