Framework for Quantifying the Efficiency of Competing Signal Transmission Modes in Proteins

TL;DR

This paper introduces a theoretical framework to quantify how efficiently different signal transmission modes work in proteins, using simulations and experiments.

Contribution

A novel theoretical framework based on linear-response theory to analyze signal transmission modes in proteins.

Findings

Signal propagation through coiled-coil motifs in histidine kinase occurs via shift, splay, and twist deformation modes.

Splay deformation is the most relevant mode for the biological function of the histidine kinase protein.

The framework was applied to analyze signal transmission in the β2-adrenergic receptor's structural domains.

Abstract

On the microscopic
level, biological signal transmission relies
on coordinated transient structural changes in allosteric proteins
that involve sensor and effector modules. The time scales and microscopic
details of signal transmission in proteins are often unclear, despite
a plethora of structural information on signaling proteins. Based
on linear-response theory, we develop the theoretical framework to
define frequency-dependent force and displacement transmit functions
through proteins and, more generally, viscoelastic media. Transmit
functions quantify the fraction of a local time-dependent perturbation
at one site, be it a deformation, a force or a combination thereof,
that survives at a distant site. They are defined in terms of equilibrium
fluctuations from simulations or experimental observations. We apply
the framework to our all-atom molecular dynamics simulations of a
bacterial histidine kinase protein extensively studied in experiments.
For the isolated coiled-coil (CC) motif that connects sensor and effector
modules, our analysis reveals that signal propagation through the
CC is possible via shift, splay, and twist deformation modes, which
is confirmed by simulations of the entire protein. Based on mutation
experiments, we infer that the most relevant mode for the biological
function of the histidine kinase is the splay deformation. For the
β2-adrenergic receptor, a transmembrane protein involved
in the G-protein signaling pathway, we compare signal transmission
across its different structural domains involved in receptor activation.

Figures

A) Picture of the full-length structure of the dark-adapted blue-light-regulated histidine kinase YF1 (PDB ID: 4GCZ), rendered using VMD. The sensor, linker, and effector modules of YF1 are homodimeric at the molecular level. (B) Schematic representation of signal transmission from the sensor to the effector site. (C) Simulation unit cell containing only the coiled-coil (CC) linker. Water and ions are not shown. (D) Time series and corresponding probability density functions (PDFs) for the N and C-termini distances of the isolated CC, |n⃗| and |c⃗|. (E) Schematic representation of shift, splay, and twist modes of the CC. (F) Schematics for defining the shift vectors n⃗ sh, c⃗ sh and the splay vectors n⃗ sp, c⃗ sp as the parallel and perpendicular components of the distance vectors n⃗, c⃗ with respect to the symmetrized CC axis a⃗ = (a⃗ 1 + a⃗ 2)/2. (G) Schematic for describing twist as rotation of splay unit vectors n̂ sp and ĉ sp around the long-axis a⃗. The angle between n̂ sp and ĉ sp is denoted as the twist angle ϕ.

(A) Self- and cross-correlation functions (C self and C cross) for the twist mode of the isolated CC. (B) Response functions, J, obtained using eq from numerical derivatives of the simulated correlation functions in panel A are shown as solid lines. The vertical dash-dotted line represents the cutoff beyond which the response functions are set to zero to prevent noise artifacts when calculating Fourier transforms. (C) Real and (D) imaginary parts of J̃(ω) obtained from the discrete Fourier transform of J(t) are shown as solid lines. Dashed lines represent simultaneous fits to the real and imaginary parts by a sum of 10 Debye relaxation functions (for details, see Section S4 in the SI). Dashed lines in panel B represent the inverse Fourier transform of the fits in C and D.

(Top) Real and (bottom) imaginary part of force transmit functions T̃ F for the shift (sh), splay (sp), and twist (tw) modes of the isolated CC obtained using analytical representations for the self- and cross-response functions (see details in Section S4 in the SI) according to eq . The sensor-to-effector (s → e) and effector-to-sensor (e → s) transmit functions are shown as solid and dashed lines, respectively.

Transmitted force profiles for the different modes of the isolated CC for different input force signals: (A) δ pulse, (B) step function, and (C) rectangular pulse of width τ = 1 ns. Insets in panels A and C represent zoomed-in force profiles. The force profiles in B and C are rescaled by the input force strength.

Effects of single-point mutations, Q133L (red), R135L (green), and charged NH3+ N-termini (blue) on the step force transmit function, F stp(t), of the isolated CC via (A) shift, (B) splay, and (C) twist modes. For comparison, F stp for each deformation mode is shown from Figure B for the wild-type CC with charge-neutral termini, denoted as AcNH–CC–CONH2 (black). Solid lines and broken lines (only shown for the wild-type CC with charge-neutral and charged N-temini) represent step responses at the effector side, Festp(t) , and at the sensor side, Fsstp(t) , respectively. All force profiles are rescaled by the input force strength.

Signal transfer in the GPCR β2-adrenergic receptor (β2-AR) protein. A) Active and B) inactive state structures of β2-AR, for the whole simulation box see Figure S9A,B in the SI. The intracellular part of the transmembrane helix 6 (TM6) remains in the open or close conformation in the active or inactive state, respectively. Signal transmission is characterized by monitoring the two extracellular separationsE1 (distance between center-of-mass of residues D113 and Y316) and E2 (distance between the S207 side-chain oxygen and one oxygen of the agonist’s catechol group)and the two intracellular separationsI1 (distance between center-of-mass of the side chains of Y219 and Y326) and I2 (distance between the backbone Cα atoms of R131 and L272). Modulus of the force transmit functions, |T̃ F(ω)|, for signaling between different extracellular and intracellular parts of β2-AR in the C) active and D) inactive state.

Tables

1. Correlation Coefficients between Different Deformation Modes for the Same Ends (Sensor and Effector) and for Different Ends (Cross) of the CC from the Simulation of the Whole Histidine Kinase Protein

Type	sh-sp	sp-tw	sh-tw
sensor	0.19 ± 0.06	0.10 ± 0.05	0.03 ± 0.09
effector	0.60 ± 0.05	-0.05 ± 0.12	0.03 ± 0.08
cross	0.31 ± 0.06	0.09 ± 0.10	0.05 ± 0.05

Full text

Introduction

Signal transmission within a cell, between cells, and from the exterior to the interior of a cell is necessary for any form of life, be it bacteria, plants, or animals. Signal-transducing units, at the molecular level, involve signal receptors (sensors) that sense intracellular or environmental changes (such as photons, hormones, or pH) and corresponding effectors that are responsible for sparking a response. ?,? Sensor and effector modules are usually distinct domains of membrane-bound or cytosolic proteins. ?−? ? How information, at the molecular level, is transmitted through proteins has been the subject of intense research in the last few decades. ?−? ? ? ? ? ? ? ? Experimental techniques, such as NMR spectroscopy ?,? time-resolved crystallography ?,? cryo-electron microscopy,? and single-molecule experiments ?−? ? have provided information about protein structure, dynamics, and mechanical signaling pathways. Computational approaches combining molecular simulation techniques ?−? ? ? ? and tools from information theory? and graph theory ?−? ? along with various linear or nonlinear correlation analyses ?,? have provided molecular-level insights into protein allosteric communication pathways. ?,? A notable mention is the statistical mechanical approach to quantify the allosteric coupling between two distal regions in a protein via the thermodynamic coupling function.? Many of the theoretical frameworks, however, use symmetric functions, which in principle are insufficient to quantify the directional, source-to-target transfer of bimolecular information? and provide only stationary (time-independent) information about allosteric coupling. A combination of spectroscopy techniques and molecular dynamics (MD) simulations has led to recent insights into time-resolved processes of allosteric regulation. ?−? ? ? ? ? However, the quantitative relation between the dynamics of a protein and its signal transfer efficacy, which ultimately governs physiological response, is missing.

We show how to quantify the signal-transfer efficiency through proteins in terms of frequency-dependent force and displacement transmit functions. A force transmit function describes the fraction of the frequency-dependent force applied at the protein sensor position that survives at the effector position. The displacement transmit function is defined similarly but is based on spatial displacements. In fact, the presence of strongly correlated fluctuations at the sensor and effector positions, as quantified by two-point correlation functions, is not sufficient for efficient signal transmission, because equilibrium fluctuations produce a background that the signal has to compete with. Transmit functions weigh the correlation between the sensor and effector positions by the fluctuation magnitude and therefore quantify the efficiency of signal transmission. They are thus distinctly different from ordinary two-point correlation functions. Our theoretical framework builds on our previously developed convolution theory? and is exact on the linear-response level.? From the time-domain transmit function, the response of a protein to any temporal perturbation signal can be calculated by convolution, all one needs as input to our theory are time series of positions or displacements that can be obtained from MD simulations or single-molecule experiments, e.g., by fluorescence resonance energy transfer.?

We perform all-atom MD simulations of an engineered blue-light-regulated histidine kinase protein.? In this synthetic enzyme, a coiled-coil (CC) motif connects the light-oxygen-voltage sensor module from Bacillus subtilis YtvA and the histidine kinase complex from the bacterium Bradyrhizobium japonicum, see FigureA. CCs are frequently found in various plant and bacterial signal transduction systems, that connect signal receptor and response or effector proteins. ?,?−? ? ? ? ? ? Signal transducers are commonly oligomeric domains and modular in architecture, produced by recombination of sensor and effector modules. Signaling through CCs leads to conformational rearrangements within effectors that trigger interaction with regulator proteins, starting signaling cascades. Studies investigating the role of CCs in signal propagation have, however, been largely limited to quantifying their stationary structural changes, such as helix pivoting or rotation ?,? in the presence of external stimuli. We study the dynamics of the isolated CC and the CC connected to the sensor and effector modules in terms of its shift, splay, and twist deformation modes. Computed transmit functions reveal that all modes are able to transmit signals. From the time-domain transmit functions, we derive transmit properties of different time-dependent signals. For a step signal, we find that transmission via the splay mode is markedly reduced by single-point mutations in the CC, whereas transmission via the twist mode is only minimally affected. Together with the experimental observation of light-induced splaying of the CC ?,? and experimental mutation studies,? this suggests that splay is the most relevant deformation mode of the CC for the biological function of this sensor histidine kinase.?

A) Picture of the full-length structure of the dark-adapted blue-light-regulated histidine kinase YF1 (PDB ID: 4GCZ), rendered using VMD. The sensor, linker, and effector modules of YF1 are homodimeric at the molecular level. (B) Schematic representation of signal transmission from the sensor to the effector site. (C) Simulation unit cell containing only the coiled-coil (CC) linker. Water and ions are not shown. (D) Time series and corresponding probability density functions (PDFs) for the N and C-termini distances of the isolated CC, |n⃗| and |c⃗|. (E) Schematic representation of shift, splay, and twist modes of the CC. (F) Schematics for defining the shift vectors n⃗ sh, c⃗ sh and the splay vectors n⃗ sp, c⃗ sp as the parallel and perpendicular components of the distance vectors n⃗, c⃗ with respect to the symmetrized CC axis a⃗ = (a⃗ 1 + a⃗ 2)/2. (G) Schematic for describing twist as rotation of splay unit vectors n̂ sp and ĉ sp around the long-axis a⃗. The angle between n̂ sp and ĉ sp is denoted as the twist angle ϕ.

To demonstrate the versatility of our method, we also apply our transmit-function formalism to a well-studied class A (rhodopsin-like) G-protein coupled receptor (GPCR), namely the β_2_-adrenergic receptor (β_2_-AR). ?−? ? ? GPCRs are the largest family of membrane proteins encoded in the human genome, responsible for detecting external signals and transducing them across the cell membrane via activation of intracellular G-proteins. ?−? ? These heptahelical transmembrane receptors are involved in virtually every physiological process and presently are targets for approximately 36% of all approved drugs in modern medicine.? β_2_-AR regulates diverse physiological functions, including smooth muscle relaxation in airways and blood vessels, cardiac contractility, metabolic processes, and immune responses, and serves as a prototypical model for understanding GPCR activation, signaling pathways, and drug interactions. Our theoretical framework enables a quantitative assessment of signal propagation across β_2_-AR’s different structural domains that act as potential key allosteric hotspots during receptor activation.

Results

Theory of Transmit Functions

The input–output relation for a general responsive system can be quantified by two transmit functions: force and displacement transmission. Consider we apply to a system’s sensor (input) position a frequency-dependent force F̃ s(ω) and to its effector (output) position a force F̃ e(ω), see FigureB. In Fourier space, the sensor and effector positions X̃ s(ω) and X̃ e(ω) change to linear order as?

$$[eqn]$$

where J̃ cross(ω) and $[eqn]$ are the frequency-dependent cross and self (sensor/effector side) linear response functions, respectively. Note that there is only a single cross-response, which is a consequence of Onsager’s celebrated reciprocal relation. However, the self-responses $[eqn]$ and $[eqn]$ are different for a general asymmetric system. It should be noted that eq is exact on the linear response level and is applicable for any two positions in a protein. The force transmit function $[eqn]$ is defined for the boundary condition of a stationary effector position X̃ e = 0, from eq we obtain

$$[eqn]$$

Similarly, setting X̃ s = 0 in eq, we obtain $[eqn]$ . Response functions are causal, i.e., there is no positional response before a force is applied, which implies J̃ cross(ω) and J̃ self(ω) have no poles in the upper half of the complex frequency plane. If furthermore J̃(ω) is nonzero in the upper half of the complex frequency plane, from eq it follows that transmit functions are also causal.

The displacement transmit function $[eqn]$ is defined as the ratio of the displacement of the effector site X̃ e divided by the displacement at the sensor site X̃ s under force-free boundary condition, F̃ e = 0. Inverting eq yields

$$[eqn]$$

where G̃’s are the moduli determined by inverting the response matrix $[eqn]$ . From eq, we obtain

$$[eqn]$$

Similarly, for F̃ s = 0, we obtain $[eqn]$ . Thus, force and inverse displacement transmit functions are the same. Note the striking resemblance between the transmit functions defined here and the transfer function which characterizes the output of a linear time-invariant system (e.g., an electric circuit consisting of resistors, inductors, and capacitors) in the context of signal processing.?

In practice, one need not apply external forces to determine linear response functions. J(t) can be obtained from equilibrium time-correlation functions, C(t), using the fluctuation–dissipation theorem? (for a derivation, see Section S1 in the Supporting Information (SI))

$$[eqn]$$

where k B is the Boltzmann constant, T represents temperature, and θ­(t) is the Heaviside step function. The needed cross and self-correlation functions are defined as C cross(t) = ⟨X s(0)X e(t)⟩ and $[eqn]$ , $[eqn]$ , respectively (see Methods for details). Note that the positional trajectories of two sites within a protein, required for the calculation of C(t), can be obtained not only from MD simulations but also from single-molecule experiments.?

Structural Stability of

the CC Linker from MD Simulations

We consider the signaling protein histidine kinase introduced above to demonstrate the applicability of our theoretical framework to real systems. We perform explicit solvent all-atom MD simulations of the whole protein for the wild-type sequence (simulation details provided in Methods, and the simulation box shown in Figure S5 in the SI) and of the isolated CC linker for the wild-type and two different mutants. We mainly concentrate on the analysis of the CC-only simulations because of the superior data quality (20 μs-long simulation each) and validate our obtained results by comparison with the whole-protein simulation. For the isolated CC systems, a summary of simulation details is given in Methods, and the simulation unit cell is shown in FigureC. We find that the distribution of the N-termini distance at the sensor side is broader compared to the C-termini distance at the effector side, both distributions exhibit tails that reflect intermittent splaying of the α-helix termini (see FigureD,E). To check the long-time stability of secondary and tertiary structures of the isolated CC, we calculate three different order parameters: the fraction of native contacts Q between the two α-helices, the root-mean-square deviation (RMSD) of distances between the native and simulated CC structures (the former taken from the crystal structure of the full-length protein shown in FigureA), and the secondary structure (SS) content. These order parameters are defined in Section S2 in the SI, and their time-averaged values are provided in Figure S1. The overall configuration of the CC remains stable within 20 μs of simulation as indicated by an average RMSD of 2.5 Å. The two α-helices remain bound to each other (Q = 0.92) due to salt bridges (involving residues R135, E138, and E142) and hydrophobic interactions (involving residues L136, L139, L143, and V146). ?,? Individual α-helices also remain stable, as their fraction of SS content values exceed 0.85.

The CC Linker Transmits Signals via Shift,

Splay, and Twist Modes

We consider CC linker deformations that result from forces of equal magnitude and opposite direction acting on the two N-termini and the two C-termini. Thus, these deformations conserve linear momentum. We define three distinct deformation modes by different orientations of the terminal displacement vector with respect to the distance vector between the terminal groups: splay, where the displacement is parallel to the terminal separation, and shift and twist, where the displacement is perpendicular to the terminal separation (see schematics in FigureE). Splay conserves angular momentum, whereas shift and twist do not and are counteracted by a rotation of the sensor module. Since the typical rotational diffusion time of the sensor module, estimated to be of the order of 1 μs (see Methods), is much longer than shift and twist relaxation times, as will be shown below, shift and twist modes are nevertheless possible signal transmission modes.

The three signal transmission modes are obtained as follows. The position vectors of the N-termini and C-termini are N⃗ i and C⃗ i, respectively, where i = 1, 2 refers to the first and second α-helix. These vectors are used to construct the separation vectors c⃗ = C⃗ 2 – C⃗ 1, n⃗ = N⃗ 2 – N⃗ 1, and a⃗ i = C⃗ i – N⃗ i (FigureF). From the end-to-end vectors a⃗ i, we define the long axis of the CC as a⃗ = (a⃗ 1 + a⃗ 2)/2. With respect to a⃗, we separate n⃗ into the parallel component, shift (n sh = |n⃗ sh| = |n⃗ · â|), and the perpendicular component, splay (n sp = |n⃗ sp| = |n⃗ – n⃗ sh|). Similarly, for the C-termini we obtain c sh and c sp. The twist angle ϕ, depicted in FigureG, is obtained from the scalar product of the two splay unit vectors n̂ sp and ĉ sp as ϕ = cos^–1^(n̂ sp · ĉ sp). We show time series of these deformation modes and the corresponding probability distribution functions in Figure S2 in the SI. Importantly, we find a strong bias toward clockwise rotation of the CC around its long axis a⃗, from the handedness plot obtained from the simulation (see Figure S2 in the SI), in agreement with the experimental observation of the left-handed supercoiling of the CC.?

The N and C-termini of the CC correspond to the sensor (s) and effector (e) side, respectively. For the shift mode, the two self-correlation functions are defined as $[eqn]$ and $[eqn]$ and the cross-correlation function is defined as $[eqn]$ . Correlation functions for the splay mode are defined by interchanging the shift quantities with the related splay quantities, e.g., $[eqn]$ . For the twist mode, the cross and two self-correlation functions are defined by the scalar product of the splay unit vectors as $[eqn]$ and $[eqn]$ and $[eqn]$ , respectively. To disentangle twist from overall CC rotation, we calculate twist correlation functions in the molecular coordinate frame obtained by removing the CC center-of-mass translation and rigid-body rotation around its principal axes at each time step.

We present results for the twist mode of the isolated CC in Figure and for shift and splay modes in Figure S3 in the SI. Self-and cross-correlation functions are shown in FigureA, the corresponding response functions J(t), obtained using eq, are shown in FigureB. All self-and cross-response functions smoothly decay to zero. The relaxation time, defined as the largest decay time τ_max_ of a multiexponential fit of J (Section S4 in the SI), is found to be the fastest for twist (τ_max_ = 2.6 ns), followed by shift (τ_max_ = 7.1 ns) and splay (τ_max_ = 9.4 ns). As expected, for each mode, the two self-response functions are greater than the cross response at all times. The real and imaginary parts of the Fourier-transformed response functions, Re J̃(ω) and Im J̃(ω), are shown in FigureC,D. There is a distinct low-frequency plateau/peak in the real/imaginary part of the twist response. To obtain analytical representations, we fit multi-Debye functions (dashed lines) to the Fourier-transformed response functions J̃(ω) in FigureC,D (details are provided in Sections S3 and S4 in the SI). Inverse Fourier transforms of the fit functions (dashed lines) reproduce the simulated time-domain response functions (solid lines), J(t), in FigureB very well.

(A) Self- and cross-correlation functions (C self and C cross) for the twist mode of the isolated CC. (B) Response functions, J, obtained using eq from numerical derivatives of the simulated correlation functions in panel A are shown as solid lines. The vertical dash-dotted line represents the cutoff beyond which the response functions are set to zero to prevent noise artifacts when calculating Fourier transforms. (C) Real and (D) imaginary parts of J̃(ω) obtained from the discrete Fourier transform of J(t) are shown as solid lines. Dashed lines represent simultaneous fits to the real and imaginary parts by a sum of 10 Debye relaxation functions (for details, see Section S4 in the SI). Dashed lines in panel B represent the inverse Fourier transform of the fits in C and D.

Force transmit functions, T̃ F(ω), obtained from the fitted J̃ self(ω) and J̃ cross(ω) using eq are shown for all three different signaling modes in Figure. T̃ F(ω) quantifies the system’s response to all possible excitation frequencies. We find that for the entire frequency range, force transmission through the twist mode (blue) is the highest, followed by the shift mode (red) and the splay mode (green). The differences are caused by energetic as well as dissipative effects, the latter include internal friction as well as hydrodynamic friction with the embedding solvent. From Re T̃ F(ω) in Figure, it is also evident that no force transmission is possible via the shift and splay mode for an input signal of frequency >0.4 THz (ps^–1^) and >0.01 THz, respectively. These cutoff frequencies are similar to the water Debye mode at a frequency of about 0.02 THz,? which suggests that the damping is partially due to the coupling to the hydration water. Except for the twist mode, the transmit functions are generally asymmetric, i.e., the sensor-to-effector side transmit function $[eqn]$ and the effector-to-sensor side transmit function $[eqn]$ are different. It should be noted that T̃ F(ω) presented in Figure characterizes the transmit properties of an isolated CC. The effects of added sensor and effector protein modules are obtained from the whole histidine kinase protein simulation (for details, see Section S6 and Figures S5–S7 in the SI). Although the absolute values of T̃ F(ω) obtained from the simulation of the whole protein are larger than those of the isolated CC, the trends for the different signaling modes remain the same. We conclude that the reduced CC terminal fluctuations in the full construct in fact increase the signal transmission efficiency compared to the isolated CC scenario.

(Top) Real and (bottom) imaginary part of force transmit functions T̃ F for the shift (sh), splay (sp), and twist (tw) modes of the isolated CC obtained using analytical representations for the self- and cross-response functions (see details in Section S4 in the SI) according to eq . The sensor-to-effector (s → e) and effector-to-sensor (e → s) transmit functions are shown as solid and dashed lines, respectively.

Different force-field parameters can significantly alter protein stability and dynamics in simulations. We check this by performing simulations of the isolated CC using Amber99SB-ILDN ?,? and DES-Amber ?−? ? in addition to CHARMM36m? used for the results presented in the main text. We find the signal transmission efficiency to be rather robust with respect to force-field variations, see Section S7 in the SI for details.

Coupling between Deformation

Modes of the CC

Berntsson et al.? have experimentally observed superhelical coiling of the CC subsequent to its light-induced splaying apart of the histidine kinase sensor domains, which suggests coupling between splay and twist deformations.? We quantify the coupling between different deformation modes by calculating the Pearson correlation coefficients (see Methods) using time-series data from the simulation of the whole histidine kinase protein, as summarized in Table. We observe that the shift-splay coupling is the highest, followed by the splay-twist coupling, regardless of whether considering data for the same ends (sensor and effector) or the different ends (cross). The shift-twist coupling is found to be insignificant in all cases. We indeed find nonvanishing splay-twist coupling for the sensor side, which explains the experimentally observed twisting by coupling to splay.

1: Correlation Coefficients between Different Deformation Modes for the Same Ends (Sensor and Effector) and for Different Ends (Cross) of the CC from the Simulation of the Whole Histidine Kinase Protein

Signal Transmission in the Time Domain

For practical purposes, signal transmission is conveniently characterized in the time domain. The force transmit function in the time domain, $[eqn]$ , is given by the inverse Fourier transform of $[eqn]$ and describes the transmission of a δ-function input force signal F s(t) at the sensor side. Once $[eqn]$ is known, the transmitted force F e(t) at the effector side due to an arbitrary input force signal F s(t) can be obtained via convolution,

$$[eqn]$$

T F(t) for the different signaling modes is shown in FigureA, and the method to obtain analytical representations for T F(t) from T̃ F(ω) is explained in Section S5 in the SI. We find that the transmitted force signals for all the different modes decay rather quickly over a few picoseconds, which reflects the dominant fast time scale of the relaxation spectrum (see Tables S2–S4 in the SI).

Transmitted force profiles for the different modes of the isolated CC for different input force signals: (A) δ pulse, (B) step function, and (C) rectangular pulse of width τ = 1 ns. Insets in panels A and C represent zoomed-in force profiles. The force profiles in B and C are rescaled by the input force strength.

In reality, signals are not transmitted via infinitely short δ-pulses but rather by pulses of finite duration. How efficiently the CC linker transmits such such a finite-duration input force signal depends on how quickly it responds to a suddenly imposed force and how quickly it relaxes back to its equilibrium state after the removal of force. To understand these complex dynamics, we discuss the transmitted output forces for force-step and rectangular-force-pulse input signals, $[eqn]$ and $[eqn]$ , both obtained from the convolution integral eq (for details, see Section S5 in the SI). The step-response function is relevant for understanding activation of signaling proteins, since an agonist binding leads to a finite, sustained increase in the volume of the binding pocket, which can be modeled as a step-like increase in the separation between any two protein residues near the binding pocket. For a force switched on at t = 0, $[eqn]$ is shown in FigureB for the different modes. The transmitted force rises steeply on the subpicosecond scale but stationary plateau values are reached for the different modes after 30–50 ns, which reflects the longer time scales of the relaxation spectrum (see Tables S2–S4 in the SI) and is around 2 orders of magnitude faster than the reported experimental time scale of 2 μs associated with light-induced conformational transitions within the sensor module.? Signal transmission through all three modes is thus sufficiently rapid to avoid becoming rate-limiting. The plateau value of F e ^stp^ is the highest for the twist mode, followed by shift and splay modes. In FigureC, we show the transmitted force $[eqn]$ for a rectangular force pulse signal of duration τ = 1 ns that is switched on at t = 0. The transmitted force through each mode decays to zero rather quickly, after the removal of the applied force. Results for different durations of the rectangular force pulse of τ = 10^–2^, 10^–1^, 10^1^, 10^2^ ns are presented in Figure S4 in the SI. Note that due to the diffusive nature of the relevant conformational transitions, which can involve rotations of entire domains, signal propagation in larger proteins can imply substantially longer time scales than found from our analysis of the CC; signal propagation time scales in general depend on the length scale and complexity of the signal-transducing system.

Robustness with Respect

to Mutations in the CC

It has been experimentally demonstrated that single-point mutations within the CC reduce the signal response of the blue-light-regulated histidine kinase YF1. ?,? To investigate this using our framework, we perform MD simulations of two different experimentally studied mutants, Q133L and R135L. We find that these single-point mutations do not affect the overall coiled-coil conformation, as observed from their different structural order parameter values compared with that of the wild-type CC, provided in Figure S1 in the SI. However, the dynamics of these two mutants are completely different from each other and from the wild-type CC, as becomes clear from their step-force responses, $[eqn]$ , shown in Figure.

Effects of single-point mutations, Q133L (red), R135L (green), and charged NH3+ N-termini (blue) on the step force transmit function, F stp(t), of the isolated CC via (A) shift, (B) splay, and (C) twist modes. For comparison, F stp for each deformation mode is shown from Figure B for the wild-type CC with charge-neutral termini, denoted as AcNH–CC–CONH2 (black). Solid lines and broken lines (only shown for the wild-type CC with charge-neutral and charged N-temini) represent step responses at the effector side, Festp(t) , and at the sensor side, Fsstp(t) , respectively. All force profiles are rescaled by the input force strength.

For the mutant Q133L (red), the plateau values of $[eqn]$ for the shift and splay modes are 16.7 and 1.8 times smaller, respectively, than that of the wild-type CC (black), as shown in FigureA,B. In contrast, the plateau value of $[eqn]$ for the twist mode is 1.8 times larger than that of the wild-type CC (FigureC). Based on the absence of signaling for Q133L in experiments,? the results in Figure suggest that the signaling in the histidine kinase is not connected to the twist mode. For the mutant R135L (green), the shift and splay mode $[eqn]$ plateau values are finite but negative. Compared to the wild-type CC, the $[eqn]$ plateau value magnitudes for the shift, splay, and twist modes are reduced by 1.6, 5.7, and 1.9 times, respectively. In experiments, the signaling activity of the mutant R135L is reduced but not absent.? In addition, structural characterizations by electron paramagnetic resonance spectroscopy and X-ray solution scattering have revealed that light induces a splaying apart of the sensor domains of the histidine kinase YF1 and hence of the N-termini of the wild-type CC. ?,? By integrating these experimental findings with our comparative analysis of mode-dependent step responses shown in Figure for the wild-type and two mutated CC variantswhich exhibit absent or reduced signalingwe conclude that the biologically relevant signaling mode in the histidine kinase is predominantly of the splay type.

Interestingly, our simulations reveal the twist mode to be the most stable mode with respect to the two single-point mutations (FigureC). Thus, it is conceivable that the CC linker might in a different biological context also function as a twist transmitter.

Asymmetric Signal Transmission through the

CC

Asymmetry could be important for the efficient information transfer from the sensor to the effector side. To look into this, we introduce the rectification factor γ as the ratio of sensor-to-effector (s → e) and effector-to-sensor (e → s) step-force transmission profile plateau values,

$$[eqn]$$

From the results for the wild-type isolated CC shown in FigureB, we conclude that rectification for the splay mode is the highest, γ^sp^ = 3.2, followed by the shift mode, γ^sh^ = 2.4. In contrast, no rectification is observed for the twist mode, i.e., γ^tw^ = 1. These results can be rationalized by noting that the rectification factor is essentially determined by the ratio of the real parts of the low-frequency self-responses for the sensor and the effector ends (Figure S3C in the SI). To study the relation between the rectification factor and the structural asymmetry in more detail, we introduce an additional asymmetry between the sensor and effector ends of the CC by uncapping the sensor side α-helix termini, which thereby become positively charged at neutral pH, resulting in the structure $[eqn]$ –CC–CONH_2_ (note that the results presented in Figures–? are obtained for the CC linker with charge-neutral end groups: AcNH–CC–CONH_2_). Though the step-force transmission profiles for AcNH–CC–CONH_2_ and $[eqn]$ –CC–CONH_2_ are qualitatively the same, we observe a pronounced difference in the plateau values of $[eqn]$ and $[eqn]$ for the shift and splay modes of $[eqn]$ –CC–CONH_2_ (see FigureA,B). However, the $[eqn]$ and $[eqn]$ plateau values for the twist mode are the same and remain unaffected in comparison to that of the charge neutral-termini system, AcNH–CC–CONH_2_ (FigureC). We thus find that the rectification factor γ can be tuned by changing chemical structures of sensor and effector-side terminal groups.

For the whole histidine kinase protein, which exhibits an additional asymmetry due to the added sensor and effector modules, we find rather similar rectification factor for the splay mode, γ^sp^ = 3.63, in comparison to the isolated wild-type CC. The corresponding value for twist (γ^tw^ = 0.82) is slightly reduced, while that for the shift (γ^sh^ = 0.67) is significantly reduced, compared to the isolated wild-type CC. Thus, the rectification factor for the splay mode is significantly larger than that for the shift and twist modes for the whole protein as well. Interestingly, this may be closely associated with the light-induced splaying of the CC linker observed in experiments. ?,? This observation suggests a possible correlation between the γ values of the different modes and their respective functional relevance.

Signal Transmission in the GPCR β2-AR Protein

We have, so far, discussed signal propagation in the histidine kinase YF1 from the sensor to effector module through the CC linker, in which the two α-helices primarily exhibit shifting, splaying, and twisting. Cytosolic or transmembrane proteins, in general, show more complex deformations including melting and unfolding of helices. Our transmit function formalism is also applicable to the study of such cell signaling proteins with complex topologies, such as GPCRs having seven helical transmembrane domains. ?−? ? We consider a prototypical GPCR, the β_2_-AR protein (FigureA,B), which plays a crucial role in regulating cardiovascular, pulmonary, and metabolic functions. As a key target for bronchodilators in the treatment of asthma and chronic obstructive pulmonary disease,? β_2_-AR is also implicated in biased signaling, where different ligands selectively activate distinct intracellular pathways, offering therapeutic potential beyond traditional agonists.? Structural insights into β_2_-AR have significantly advanced our understanding of GPCR activation, ligand specificity, and allosteric modulation.? Recent computational and cryo-electron microscopy studies have provided insights into the activation and deactivation mechanisms of β_2_-AR. ?−? ?,?

Signal transfer in the GPCR β2-adrenergic receptor (β2-AR) protein. A) Active and B) inactive state structures of β2-AR, for the whole simulation box see Figure S9A,B in the SI. The intracellular part of the transmembrane helix 6 (TM6) remains in the open or close conformation in the active or inactive state, respectively. Signal transmission is characterized by monitoring the two extracellular separationsE1 (distance between center-of-mass of residues D113 and Y316) and E2 (distance between the S207 side-chain oxygen and one oxygen of the agonist’s catechol group)and the two intracellular separationsI1 (distance between center-of-mass of the side chains of Y219 and Y326) and I2 (distance between the backbone Cα atoms of R131 and L272). Modulus of the force transmit functions, |T̃ F(ω)|, for signaling between different extracellular and intracellular parts of β2-AR in the C) active and D) inactive state.

We analyze two different simulation trajectories (see Methods) of a high-affinity agonist-bound β_2_-AR protein, embedded in a lipid membrane, in the (i) active and (ii) inactive state. ?,? The β_2_-AR transitions from its active to inactive state following decoupling of the G protein, accompanied by structural rearrangements that stabilize the inactive conformation.? To evaluate whether our transmit function formalism can pinpoint critical hotspots of β_2_-AR’s allosteric signaling pathways, we focus on reaction coordinates known to play important roles in the receptor activation process. ?−? ? Those are the two extracellular separations, E1 and E2, and the two intracellular separations, I1 and I2, as depicted in FigureA,B. The selection of the reference atoms has been guided by functional considerations in addition to structural stability, as described below. E1 (D113–Y316 distance) and E2 (S207–agonist distance) monitor conformational changes in the ligand binding site. As a matter of fact, S207 is a key residue in the orthosteric binding pocket of β_2_-AR, where it forms direct hydrogen bonds with catechol ligands such as epinephrine, and plays a central role in receptor activation.? I1 (Y219–Y326) is an important contact pair since reorientation of the side chains is functionally relevant and provides a meaningful indicator of gate opening on the intracellular side, in addition, Y326 is conserved in 99% of class A GPCRs. I2 (R131–L272 distance) captures the outward movement of transmembrane helix 6 (TM6) and the opening of the intracellular gate. Specifically, R131 belongs to the highly conserved DRY motif that stabilizes the inactive state and regulates conformational changes during activation. ?,? L272, located in TM6, has been deliberately selected as a reference point because the outward movement of TM6 at the cytoplasmic site is the hallmark of GPCR activation and is critical for G-protein coupling.?

Regardless of the different protein conformational ensembles sampled in the (i) active and (ii) inactive states, our transmit-function analysis shows that signaling from E2 to I2 is the most efficient, as is clear form the low-frequency plateau values of the modulus of the force transmit functions |T̃ F(ω)| shown in FigureC,D. The efficient signaling between E2 (S207–agonist) and I2 (R131–L272) is consistent with a well-supported scenario based on experimental findings, according to which ligand engagement at S207 initiates conformational changes that propagate through the DRY motif to the cytoplasmic opening, where TM6 displacement exposes the intracellular space for G-protein binding. ?,?−? ? We thus demonstrate that our transmit-function formalism allows us to compare the efficiency of different signaling pathways in complex systems and to find optimal pairs of sensor and effector sites.

Discussion and Conclusions

We present the theoretical framework to quantify the signal transmission between two distinct sites of a protein through different deformation modes, expressed in terms of the corresponding self- and cross-response functions. According to the fluctuation–dissipation theorem, the response functions are related to equilibrium time-correlation functions, which can be obtained from MD simulations as well as from single-molecule experiments. Note that in experiments, trajectories of separation coordinates typically include effects due to the coupling to measurement devices, which can be removed by using dynamic deconvolution theory. ?,? The displacement transmit function relates the correlations between two sites and the fluctuations at the site at which the input signal is applied, it thus quantifies the ratio of the output to the input signal and thus conveys more useful information than the often considered dynamic cross-correlation.

Applying our theoretical framework to the histidine kinase protein YF1, we demonstrate that all three deformation modestwist, shift, and splayin principle enable signal transfer from the sensor to the effector end of the CC. The experimentally observed splaying followed by superhelical coiling of the CC upon activation? is expected due to the coupling between splay and twist deformation modes, as quantified by our analysis. Analysis of our simulation data for the wild-type CC and two single-point mutants? suggests that splay is actually the signaling mode realized in the experimentally studied histidine kinase. ?,? Although twist, in principle, is a more effective mode of signal transmission, it does not conserve angular momentum and would therefore induce rotation of the sensor domain, this is probably why nature is not using it, at least for this protein construct. Previous experiments have indicated that the length of the CC linker, not only the actual linker sequence, is instrumental in determining the response to light signals.? Our framework could in the future also be used to study the CC length-dependent signaling. In this paper, we have focused on signaling between identical deformation modes at the CC termini: twist-to-twist, splay-to-splay, and shift-to-shift. Off-diagonal signaling modes, that means coupling of e.g., splay at the sensor side to shift at the effector side, could be experimentally relevant and will be considered in future work.

For the GPCR β_2_-AR protein, our transmit-function formalism enables quantitative characterization of signal transmission across its structural domains, pinpointing critical allosteric sites implicated in activation. Our framework, with its capability to quantitatively map signal propagation between different domains, offers a versatile approach to study allosteric signal propagation in GPCRs, where complex topologies and conformational transitions challenge conventional analyses.? In future work, it could be applied to quantify how ligand binding, G-protein coupling and decoupling, or recruitment of arrestin, which are GPCR-regulating proteins, modulate signal transmission across the receptor’s structural domains.

Signal transmission through general protein networks can be predicted from the response functions of individual components by repeated application of convolution relations for serial and parallel connections. ?,? Our study, thus, provides a way forward to relate atomic-level protein dynamics to large-scale intermolecular communications of biological signaling networks.

Our theory is formulated at the linear-response level and is therefore scale-invariant with respect to the input signal amplitude. To obtain the signal threshold beyond which the signal strength surpasses the noise background, one needs to compare the signal strength with the root-mean-square of the fluctuating force or displacement, similar to the definition of the signal-to-noise ratio in information theory.?

Methods

Models and

Force-Field Parameters

Isolated CC Systems

From the crystal structure of the histidine kinase protein YF1 (PDB ID: 4GCZ), sensor (N-terminal) and effector (C-terminal) modules are deleted to obtain the structure of the isolated CC linker.? The CC is composed of two parallel α-helices, each containing the same 23 residues ([126]­Ile-Thr-Glu-His-Gln-Gln-Thr-Gln-Ala-Arg-Leu-Gln-Glu-Leu-Gln-Ser-Glu-Leu-Val-His-Val-Ser-Arg[148]). The CC is simulated in a rhombic dodecahedron box of volume 227 nm^3^ filled with 7135 water molecules (and counterions needed to neutralize the system). CHARMM36m protein force field parameters,? the TIP3P water model ?,? and ion parameters from ref.? are used. Four different systems with changes of N-termini capping groups or a mutated residue are considered: $[eqn]$ –CC–CONH_2_, AcNH–CC–CONH_2_, AcNH–CC­[Q133L]–CONH_2_, AcNH–CC­[R135L]–CONH_2_. The acetyl (AcNH) group is used at the N-terminal and the “–CONH_2_” group is used at the C-terminal Arg, to simulate charge-neutral termini. The two mutated systems are selected from the study by Gleichmann et al. ? For AcNH–CC–CONH_2_, two additional simulations using different protein force fields, Amber99SB-ILDN ?,? and DES-Amber ?−? ? are performed (for details, see Section S7 in the SI).

The Whole

Histidine Kinase Protein

The histidine kinase YF1 (PDB ID: 4GCZ), with two flavin mononucleotide (FMN) cofactors bound to its sensor module and one adenosine diphosphate (ADP) bound to its effector module, is used as the initial structure. Missing hydrogen atoms are added using CHARMM-GUI.? FMN is parametrized using the CHARMM General Force Field (CGenFF),? while the protein and ADP are assigned parameters from the CHARMM36m force field.? The complex is placed in a 15 × 15 × 15 nm^3^ cubic box and solvated with 110,941 TIP3P water molecules ?,? as shown in Figure S5 in the SI. The system is charge neutralized, and 0.10 M NaCl? is added using the Monte Carlo ion placement method in CHARMM-GUI,? ensuring a physiologically relevant ionic strength.

The GPCR β2-AR Protein

Details of the model building and equilibrium MD simulation protocol for the all-atom explicit solvent simulations of a high-affinity agonist, conformationally constrained epinephrine (c-Epi)-bound β_2_-AR protein, embedded into a pure 1-palmitoyl-2-oleyl-sn-glycero-3-phosphocholine (POPC) bilayer membrane, coupled to a GDP-bound G_s_ protein can be found in the recent publications by one of us. ?,? The CHARMM36m? parameters for protein, the CHARMM36? parameters for lipid, the TIP3P water model ?,? the parameters for c-Epi and GDP derived from CGenFF,? and ion parameters from ref.? are used. For our analysis, we consider from the above studies two different simulation trajectories of the agonist-bound β_2_-AR in the (i) active and (ii) inactive state, taken from simulations with the initial conformations denoted as structure 16 and structure 20, respectively, according to Papasergi-Scott et al.? The simulation box for each of the two states is depicted in Figure S9 in the SI. Each simulation trajectory is of adequate duration (3.15 μs) with the data saved every 10 ps. For the simulation started with structure 20 that leads to the inactive state as monitored by the TM6 distance (I2) shown in Figure S9C in the SI, the trajectory excluding the first 200 ns, which corresponds to the transition from active-to-inactive state, is used for analysis.

MD Simulation

Details

For the whole histidine kinase protein and each isolated CC system, simulations are performed in the NpT ensemble at temperature T = 300 K and pressure p = 1 bar with periodic boundary conditions using the Gromacs package? for 2 and 20 μs, respectively. The stochastic velocity rescaling thermostat? with a time constant τ_ T _ = 0.1 ps is used to control temperature, while for pressure control an isotropic Parrinello–Rahman barostat? is used with a time constant τ_ p _ = 2 ps and compressibility κ = 4.5 × 10^–5^ bar^–1^. The LINCS algorithm? is used to constrain the bonds involving hydrogen atoms, allowing a time step Δt = 2 fs. Electrostatic interactions are computed using the particle mesh Ewald method? with a real-space cutoff distance of 1.2 nm, while van der Waals interactions are modeled using Lennard-Jones potentials with a cutoff distance of 1.2 nm where the resulting forces smoothly switch to zero between 1 to 1.2 nm. For data analysis, simulation trajectories are saved every 1 ps. Images are rendered using the visual molecular dynamics (VMD) software.? Analysis is performed using in-house developed codes and Gromacs analysis modules.? The time-correlation function for observables A(t) and B(t) is calculated as

$$[eqn]$$

where τ is the time lag, and L is the trajectory length.

Protein

Rotational Relaxation Time

The rotational relaxation time τ_r_ of an object is related to the rotational diffusion coefficient as D r = 1/2τ_r_ and is estimated from the Stokes’ rotational diffusion coefficient $[eqn]$ .? Using the viscosity of the medium as that of water, η = 8.9 × 10^–1^ Pa·s, and the hydrodynamic radius as half of the largest length scale of the full-length protein YF1, R h = 7 nm, we obtain τ_r_ = 2 μs, which is notably longer than the transmission relaxation times, by 1–2 orders of magnitude.

Coupling between Deformation

Modes

Pearson correlation coefficients between two different deformation modes, i and j, are obtained as

$$[eqn]$$

where ⟨·⟩ denotes the time average, and x _ i/j _ refers to the time-series data for shift, splay, or twist at the same ends or at different ends of the CC. Using the time series data obtained from the simulation of the whole protein YF1, the computed correlation coefficients R _ ij _ for i ≠ j (as R _ ii _ = 1) are summarized in Table.

Supplementary Material