Any periodic sequence composed of invariantly coupled degrees of freedom universally leads to dispersive wave transport. Like any other folded polymer system, a protein can be first approximated as a linear chain of AA bonded between the C _α atoms (backbone) as shown on Figure 1. However, these covalent interactions are completed by noncovalent bonds between the side chains: proteins fold, and consequently, each degree of freedom of the BB has cohesive energy that is no longer spatially homogeneous Figure 1. Based on these physical considerations, we express the lagrangian of the system as, L = ∑ i = 1 N 1 2 X ̇ i 2 − ∑ 1 ≤ i < j ≤ N 1 2 α i j ( X i − X j ) 2 (1) which allows to extract the harmonic equations of motion of such a dynamical system as ω 2 X i = ∑ j = 1 N α i j X i − ∑ j = 1 N α i j X j (2) where ω stands for the harmonic pulsation, X _i is the displacement of atom i from its equilibrium position x _i . The α _ij are constants characterizing the strength of the bond between atoms i and j. We see from this eq. that disorder (irregularity) is reflected by the quantities α _ij which account for the nature of the interactions (ex. angular or distant dependent) as well as the number of “neighbors” included in the sum. We introduce the matrix (operator): − Δ : ( δ i j ) = ∑ k = 1 N α i j , if i = j . − α i j , otherwise . (3) such that the equations of motion are casted into an eigenproblem in the form of a Laplace equation: − Δ X = ω 2 X (4)

We consider vibrations defined at the zone-edge, by the wavenumber k = π/a: X ̃ j = e i π a x j X j = ( − 1 ) j X j . We introduce a real constant C, taken slightly greater than the greatest eigenvalue ω max 2 of − Δ and write ( C − ω 2 ) X ̃ i = C − ∑ j = 1 N α i j X ̃ i + ∑ j = 1 N α i j ( − 1 ) j − i X ̃ j (5) with an operator L _h defined by: L h X ̃ = ( C − ω 2 ) X ̃ (6)

A first frequency localization landscape (LL) u _h is obtained from L h u h = ( 1 ) , (7) where (1) = (1,‥,1) ^T , which is termed high frequency LL for its ability to capture the localization of the highest frequency modes (Chalopin et al., 2019). In the case of a continuous elliptic differential operator H on a domain Ω, it has been shown that the localization landscape u defined as Hu = 1 on Ω satisfies | Ψ ( x ) | ‖ Ψ ‖ ∞ ≤ E u ( x ) where Ψ is an eigenfunction of H with eigenvalue E. We would then expect in our case to have X ̃ i ‖ X ̃ ‖ ∞ ≤ ( C − ω 2 ) u i for every 1 ≤ i ≤ N. Unlike in the continuous case this is not generally true. The maximum principle is indeed not true for the L _h operator. Components of L h − 1 are then not strictly positive and the inequality is false in general. However we numerically find that this inequality is fulfilled on most cases. Therefore, high-frequency modes are expected to be strongly localized on the maxima of u _h or equivalently on the minima of W _h = 1/u _h .

In a discrete disordered system the dispersion relation is no longer monotonous: wave-vectors at the zone edge are no longer systematically associated with the highest frequencies. For this reason low frequency modes are expected to be localized too. We introduce a complementary localization landscape that captures low frequency localization at the zone edge.

We reconsider the dynamics of the system written as: ( ω 2 + ϵ ) X ̃ i = ϵ + ∑ j = 1 N α i j X ̃ i − ∑ j = 1 N α i j ( − 1 ) j − i X ̃ j (8) and introduce the operator L _l : L l X ̃ = ( ω 2 + ϵ ) X ̃ (9) which is positive definite with ϵ a small constant. The low frequency localization landscape u _l is straightforwardly obtained from: L l u l = ( 1 ) (10)

As previously the inequality X ̃ i ‖ X ̃ ‖ ∞ ≤ ( ω 2 + ϵ ) u i is not generally true but happens to be true on most cases. On Figure 2 the effective confining potential W _l = 1/u _l (red) of the Spike protein of SARS-Cov2 is displayed. One can clearly see that low frequency modes are localized on the minima of W _l . As it can be seen on Figure 3, the 2 localization landscapes happen to be complementary: the minima of the high frequency LL are located on the same sites as the maxima of the low frequency LL and conversely. This can be explained by a simple relation between operators L _h and L _l : L h = ( C + ϵ ) I N − L l (11) where I N denotes the N × N identity matrix. L _l being symmetric and positive definite, there exists an orthogonal matrix O and a diagonal matrix Σ with strictly positive eigenvalues such that: L l = O Σ O T (12)

Note that matrix O can be chosen such that it is the matrix of change of basis from the canonical basis to the basis of envelopes at k = π/a (which are eigenvectors of L _l ). Localization landscapes are defined as: L h u h = ( 1 ) L l u l = ( 1 ) (13)

We now use the diagonalization of L _l in these two expressions: ( ( C + ϵ ) I N − O Σ O T ) u h = ( 1 ) O Σ O T u l = ( 1 ) (14) ⇒ ( O T ( C + ϵ ) I N − Σ O T ) u h = O T ( 1 ) Σ O T u l = O T ( 1 ) (15) ⇒ ( ( C + ϵ ) I N − Σ ) u h * = ( 1 ) * Σ u l * = ( 1 ) * (16) where we have denoted v* = O ^T v for any column vector v and we have used the fact that O ^T and ( C + ϵ ) I N commute. Thus we have: ( ( C + ϵ ) I N − Σ ) u h * = Σ u l * (17)

Matrices Σ and ( C + ϵ ) I N − Σ being diagonals we have for every 1 ≤ i ≤ N: ( u h * ) i = Σ i i C + ϵ − Σ i i ( u l * ) i where we have denoted (v) _i the ith component of a column vector v. The Σ _ii are eigenvalues of matrix L _l so we have Σ i i = ω i 2 + ϵ for every 1 ≤ i ≤ N. Eventually we obtain a very simple relation between components of the localization landscapes in the basis of envelopes at k = π/a: ( u h * ) i = ω i 2 + ϵ C − ω i 2 ( u l * ) i (18)

Recall that we choose constant C slightly greater than the greatest eigenvalue ω max 2 of − Δ. Therefore, for low frequency ω _i Eq. 18 implies ( u h * ) i ≪ ( u l * ) i . For high frequency ω _i Eq. 18 implies ( u l * ) i ≪ ( u h * ) i .

Here we describe how the inverse quantities for both landscapes (1/u) shall be regarded as effective confining potentials (CP). This important quantity (CP) has been priorly introduced to succesfully describe confinement of electronic wavefunctions (Arnold et al., 2016).

We demonstrate that in discrete phonon systems, any localization landscape aforementioned (e.g u _l ) also leads to the apparition of an effective potential (W _l = 1/u _l ) in the equations of motion (Eq. 25), which take a form similar to the Schrödinger equation describing a free particle (e.g electron) in an irregular quantum well.

We start with the equation defining the operator L _l : L l X ̃ = ( ω 2 + ϵ ) X ̃ , (19) which can be expanded in the continuous limit as − d i v ( A ∇ ) + V X ̃ = ( ω 2 + ϵ ) X ̃ (20)

We introduce X ̃ = u ϕ to account for the modulation of the mode with an envelop defined by u = u _l . This provides: − d i v ( A u ∇ ϕ ) − d i v ( A ϕ ∇ u ) + V u ϕ = ( ω 2 + ϵ ) u ϕ . (21)

Using the Leibniz Formula d i v ( a x ⃗ ) = a ∇ ⋅ x ⃗ + x ⃗ ⋅ ∇ a (22) we obtain − u ∇ ⋅ ( A ∇ ϕ ) − ( A ∇ ϕ ) ⋅ ∇ u − ϕ ∇ ⋅ ( A ∇ u ) − A ∇ u ⋅ ∇ ϕ + V u ϕ = ( ω 2 + ϵ ) u ϕ − u ∇ ⋅ ( A ∇ ϕ ) − 2 A ∇ u ⋅ ∇ ϕ + ϕ − ∇ ⋅ ( A ∇ u ) + V u = ( ω 2 + ϵ ) u ϕ − ∇ ⋅ ( A ∇ ϕ ) − 2 A ∇ u u ⋅ ∇ ϕ − ϕ u = ( ω 2 + ϵ ) ϕ (23) where we used the fact that − ∇⋅ (A∇u) + Vu = 1. Eventually we get − ∇ ⋅ ( u 2 A ∇ ) u 2 + 1 u ϕ = ( ω 2 + ϵ ) ϕ (24) replacing by the vector X, we obtain − 1 u l 2 div ( u l 2 A ∇ X ) + X u l = ω n 2 X . (25)

To illustrate this feature, the high frequency confinement potential has been color-coded in 3D on Figure 2.

This analogy leads to introduce the concept of effective barrier through the quantity W _h − (C − ω ²). Consequently, phonons wavenumbers become purely imaginary where (C − ω ²) − W _h < 0, that is wave amplitudes vanish with exponential decays (Figure 3 and Figure 4). This later aspect is crucial for understanding how a protein operates from the knowledge of its energy partition.

In this work, most the discussion is exemplified on the SPIKE (S) protein of SARS-CoV-2, which is a molecule that plays a key role in the virus replication (Turoňová et al., 2020). Without providing too much details, one can simply recall that it serves as a receptor to bind to the angiotensin 2 converting enzyme of the host receptors, while initiating a binding with the cell membrane before fusing with the latter (S) is formed by three chains (a trimer) of polypeptide polymers (Figure 1B). Each chain (monomer) consists of a sequence of amino acids regularly spaced by 0.38 nm (see Figure 1) while the folded structure (3D) form a very irregular object. Its dynamical properties (vibrations) play a crucial role in understanding how the virus interacts with its host (Di Paola et al., 2020). The normal modes are calculated from − Δ considering a parameter free elastic network model (pfENM-GNM) (Yang et al., 2009) as shown in Figure 3. The effective high-frequency localization potential obtained from Eq. 7 is also embedded in the figure to highlight confinement effects by the potential: just like Anderson localization, phonon confinement manifests itself by an exponential modal decay upon the effect of an effective barrier. However, W _h does not completely describe the whole spectrum (i.e., low-frequency motions). When the effective energy C − ω ² rises above the maxima of the potential W _h , novel localizing regions appear where the dynamics were priorly forbidden. Interestingly, each site that is “frozen” within a given frequency band switches to a dynamically active domain, once driven in the complementary frequency band. Figure 3 and Figure 4 illustrate this complementary relationship between W _l and W _h as predicted by Eq. 11. These potentials are said reciprocal, and the dynamical result of their interplay is termed bilocalization.

This observation leads to the next question: estimating the density of states (DOS) without solving the eigenvalue problem. Weyl’s law offers an asymptotic estimate of the DOS, and it has been successfully employed in the case of electronic wavefunctions (Arnold et al., 2016) where the number of modes with energy lower than E is estimated as N ( E ) ≈ N W ( E ) = ( 2 π ) − n ∫ k 2 + W ( x ) ≤ E ⅆ n x ⅆ n k (26) Where W corresponds to the localization landscape for electrons in a random potential and E the effective energy. In our case, since normal modes are localized in space, by Fourier properties their width in k-space is expected to be large. Thus, one may assume that all wave vectors are accessible to a state localized on a scale of few DOF. Instead of summing over the volume of phase space such that k 2 + W h ( x ) ≤ C − ω 0 2 we sum over all wave vectors of the Brillouin zone and only impose a condition on position: W h ( x ) ≤ C − ω 0 2 . To go from a continuous integral to a discrete sum we use ⅆ k → π N a , ⅆx → a and ∫ → ∑. The estimation is then: N ( C − ω 0 2 ) ≈ N W h ( C − ω 0 2 ) = 1 2 π ∑ x i , k n W h ( x i ) ≤ C − ω 0 2 π N a a = 1 2 N ∑ x i W h ( x i ) ≤ C − ω 0 2 ∑ k n 1 = 1 2 N ∑ x i W h ( x i ) ≤ C − ω 0 2 2 N = ∑ x i W h ( x i ) ≤ C − ω 0 2 1 (27)

N(x) corresponds to the number of modes having a frequency such that C − ω ² ≤ x and should be termed the eigenvalue counting function (Arnold et al., 2016). Figure 5 shows comparison of this estimation and the true counting function N ( C − ω 0 2 ) for (S).

Eq. 27 tends to overestimate the number of high-frequency modes but gives a fair approximation of the density of states. What emerges from these results, which can be generalized for any discrete disordered dynamical system, is that.

1–Disordering a molecular crystal produces two localized energy bands predicted by a double wave confinement potential. These potentials have the particularity of being complementary in the sense that (Eq. 1) temporally, they predict the locations of slow (u _l ) and fast (u _h ) motions, moreover (2) the domains associated with each temporally distinct landscapes are spatially complementary: their association reproduces the molecular domain as a whole (Figure 4).

2–The localization potential of slow motions is a reading of the local rigidity (see Sup. Mat.). Thus, any mode of effective energy higher than the quantity ∑α _ij is no longer influenced by W _l . Above a threshold corresponding to the maximum of the W _l , the localization potential W _h takes over.

3–The localized modes are subject to uncertainty Δx ⋅Δk ≥ 2π, so any strong localization of the energy of a mode implies a broadening of the spectrum associated with the spatial frequencies. Disorder effect on dispersion can be approximated by the modes’ projection at the edge of the BZ. Figure 5 illustrates this aspect.

4–Complementarity implies that whatever the characteristic time considered, for each DOF, there is always a cutoff frequency below which the dynamics will be inhibited and vice-versa. The concept of bi-localization implies that there are no frozen domains. Being active or at rest, this only occurs on a limited frequency band. This description is consistent with the equipartition theorem of energy.

Bilocalization, as it is presented here, is a theory that allows to describe in a condensed and quite robust way the dynamics of proteins, that is the evolution of the positions of the degrees of freedom in space and in time, at the two limits of the slowest and fastest time scales. This tool can be adapted to any cohesion model, atomistic or coarse-grained. If it is well known that a protein is a dynamical object with multiscale behaviors, there are in fact no tools that surpass the traditional normal mode analysis to highlight the domains characterizing the motions at a given frequency/time period. The normal mode approach widely developped in the litterature, is neither rigorous nor sufficient to establish a complete picture of proteins and enzymes dynamics. Indeed, computing eigenmodes consists in accessing to the vector space that allows to decompose each motion with a well defined orthonormal basis. In other words, the motion of each atom/AA is a linear combination of the eigenmodes. If it is possible to access to this basis by solving an eigenvalue problem, we learn nothing about the combination of these modes (i.e., we do not know the coefficients that define the linear combination, nor the phase shifts between these modes), thus the information we have access to is partial, and one cannot rigorously predict an observed motion on the sole knowledge of the eigenvectors. In addition, if we want to study large systems like the SPIKE protein or RNA-polymerase (roughly 5000 AA), it is almost unmanageable to manipulate 15 0,000 eigenvalues with eigenvectors distributed over 50,000 amino acids. In most works, the analysis is often restricted to the few lowest energy mode, which form a too small fraction of the full eigenmodes spectrum. By solving the linear system LU = 1, with the right localization operator L, i.e,. by extracting a single quantity, associated to each degree of freedom, one gets a more accurate information roughly 100 times faster to calculate because solving a linear system is much easier in terms of algoritmic complexity than calculating the full spectrum of an operator.

The analysis proposed here is not only interesting to predict rigorously the predominant motions observed in these systems, but it also allows to understand the way energy flows in the latter. That is, which domains are energetically coupled and which are isolated (Chalopin et al., 2019).

In addition, bilocalization also offers a deep description of the underlying physical mechanisms that allow to universalize the formation of a partitioning/confinement of vibrations in proteins’ domains. To this end, we have established an analogy with quantum mechanics by showing that the structural disorder of the protein, which translates into an inhomogeneity of its rigidity, produces a confinement potential for the phonons (vibration modes) of the system just as a crystal produces a potential barrier for an electron (Filoche et al.). We will see later on, how this physical description allows to unify the set of well known methods such as the analysis of proteins by contact network, the long range transmission mechanisms induced by rigidity, which are in fact metrics that have some efficiency, but do not share any common physical picture of the exchange processes involved.