Collective behavior of oscillating electric dipoles

We investigate the dynamics of a population of identical biomolecules mimicked as electric dipoles with random orientations and positions in space and oscillating with their intrinsic frequencies. The biomolecules, beyond being coupled among themselves via the dipolar interaction, are also driven by a common external energy supply. A collective mode emerges by decreasing the average distance among the molecules as testified by the emergence of a clear peak in the power spectrum of the total dipole moment. This is due to a coherent vibration of the most part of the molecules at a frequency definitely larger than their own frequencies corresponding to a partial cluster synchronization of the biomolecules. These results can be verified experimentally via spectroscopic investigations of the strength of the intermolecular electrodynamic interactions, thus being able to test the possible biological relevance of the observed macroscopic mode.

The present work is motivated by a (practical) problem of (potentially) relevant impact: the experimental confirmation or refutation of the possibility of detecting long-range electrodynamic attractive forces among biomolecules, if any. The ultimate reason for searching these electrodynamic interactions stems from the observation of the high efficiency displayed by biomolecules when moving toward their specific targets and sites of action in living cells. Biochemical players "need to know" where to go and when, and are capable to reach their cognate partners so quickly that it hardly seems to be the result of a random search driven by thermal fluctuations (Brownian motion) alone. A longstanding hypothesis surmises that in order to accelerate these encounters, selective forces acting at a long distance are needed besides standard short-range ones (like covalent bonds, van der Waals forces etc.). This mechanism of molecular recruitment at a distance could be of high relevance to biology. Unfortunately, because of technological limitations, an experimental proof or refutation of this possibility has been for a long time and is still sorely lacking. The present day technological advances allow to cope with experimental challenges that were very hard to tackle in the past. This is the case of modern methods in Fluorescence Fluctuation Spectroscopy 1,2 that we invoked in our previous studies [3][4][5] , and of Terahertz spectroscopy 6 that we suggest for the present study as a reliable experimental setup to detect collective vibrating modes emerging in identical molecules ensembles.
In particular, in order to clarify the possible activation of electrodynamic interactions, it has been previously studied how the diffusion behavior of biomolecules in solution changes depending on their concentration, via the employment of Fluorescence Fluctuation Spectroscopy [3][4][5] . Varying the concentration corresponds to varying the average intermolecular distance as a consequence of the action of surmised electrodynamic intermolecular interactions. Moreover, is it shown in 7 that, by shining a laser light on an aqueous solution of proteins a strong mesoscopic dipolar vibration for each molecule can be excited and detected in a sub-Terahertz frequency range. The detection in this case is possible thanks to dye molecules attached to each protein that harvest the incoming laser light. In particular these strong dipolar oscillations can switch-on intermolecular electrodynamic interactions possibly acting at a large distance (even up to some thousands of Angstroms).
In the following we put forward an alternative/complementary experimental strategy. To this aim, we study the dynamical behaviour emerging in an ensemble of identical molecules, each one vibrating with a low frequency mode, kept active by an external energy supply. A (mesoscopic) vibrational mode of a single molecule is the coherent vibration at the same frequency of a relevant fraction (or even of all) of its atoms, which brings about an oscillating electric dipole moment at the same frequency. The latter, in turn, activates an electrodynamic SCieNtifiC REPORTS | (2018) 8:15748 | DOI: 10.1038/s41598-018-33990-y intermolecular force whose potential decays as −1/r 3 , where r is the intermolecular distance. Then, by changing the average distance among the molecules (as a control parameter), we have observed the emergence of a collective behaviour, this time involving the whole ensemble of molecules, which is detected through neat and substantial changes of a collective observable. This observable is a collective (macroscopic) one since it depends on the dynamical behaviour of all the molecules in the ensemble, and, more precisely, it consists in the power spectrum of the time variation of a quantity proportional to the total electric dipole moment of the system. In practice, this suggests a new experimental strategy to ascertain whether electrodynamic interactions among biomolecules can be activated: if we consider aqueous solutions of biomolecules, prepared with a salted solution sufficient to shield electrostatic interactions down to a few Angstroms, then we can excite these biomolecules with a suitable external forcing (like the laser light) in order to induce a mesoscopic vibration in the system at the level of the single molecule. Then, by changing the concentration of these solutions, that is by changing the intermolecular interaction strength which is a function of the distance between the molecules, and by switching on and off a laser light shining on the solution, some significative variation could be looked for by means of spectroscopic techniques to reveal whether or not the mentioned electrodynamic interactions are present. The paper is organized as follows: in Section II the model is defined and discussed, while in Sec. III we report the outcomes of the Molecular Dynamics simulations of the model and we comment on the observed phenomenology. Section IV is devoted to some concluding remarks about the results presented throughout the present paper.

The Model
Model for the biomolecule. As already stated in the Introduction, the present work aims at understanding whether through spectroscopic experiments, presumably in the Terahertz frequency domain, an experimental confirmation or refutation can be obtained of the theoretical prediction firstly stated in 7 : whether or not it is possible to activate electrodynamic forces between biomolecules, vibrating out-of-thermal equilibrium, in aqueous solution.
In what follows, we consider a simple model for an ensemble of biomolecules randomly oriented and randomly distributed in 3D space, coupled through an interaction potential decreasing as −1/r 3 as a function of the interparticle distance r. Each biomolecule is modeled as an oscillating electric dipole composed of two material points, each of them with a mass m and the same absolute value Ze of the electric charge but with opposite sign. The positions of the positive and negative charged particles of the i-th biomolecule are respectively r +,i and r −,i . The position of the center of mass of the i-th biomolecule is indicated by R i while the direction of each dipole is Both have been considered to be fixed, so that the charged particle of each biomolecule are constrained to oscillate along their joining line, i.e. dR i /dt = 0 and = dt r d / 0 i . These assumptions may seem to be quite strong with respect to a real biological molecular system where particles both diffuse and rotate due to the collisions with the surrounding water molecules. Indeed these assumptions are justified by the fact that the characteristic time scales of the inner oscillation of biomolecules are much shorter compared with those associated with the diffusion of the biomolecules centers of mass at their rotational diffusion (See Supplementary Information, Sec. I, for a more detailed discussion). It follows that the only dynamical variables are the mutual distances between the two centers of charge of each biomolecule. The electric dipole moment is given by =t Zer t p r ( ) ( ) i i i . Despite its simplicity, this model is suited to explore the presence of collective effects on the dynamics of coupled oscillating dipoles with fixed distance and orientations representing a system of oscillating biomolecules in mutual interaction through long-range quasi-static electrodynamic field generated by their oscillatory electric dipole.
Mechanical properties of a biomolecule. The mechanical properties of each biomolecule are described by an effective potential that is supposed to exist between material charged points. A minimum of the effective potential is assumed to be attained for r i = r i0 , i.e. dV/dr i (r i0 ) = 0 and , so that the effective potential is assumed to be where the parameter Λ is the characteristic length of the oscillation amplitude for the emergence of non-harmonic contributions. The non-harmonic contribution has been included for two main reasons: firstly, it accounts for the exchange of energy of the main collective mode with other vibrational normal modes of the biomolecule; secondly, it prevents instability of the oscillations when the electric dipoles are strongly coupled among them.
Mutual quasi-electrostatic interactions among biomolecules. The physical system that we are modeling is an ensemble of oscillating biomolecules in aqueous solutions in presence of freely moving ions. Since this work aims at studying collective phenomena originated by long-range electrodynamic interactions among biomolecules, we neglect any electrostatic interaction. This assumption is well justified in presence of Debye screening which, inside living cells, has a length scale of a few Angstroms. It follows that, for the range of average intermolecular distances which is of interest here (that is, -10 10 Å 2 3 ), the contribution of electrostatic fields is negligible. To the contrary, electrodynamic fields of sufficiently high frequency are not screened in water also in presence of freely moving ions, as it follows both from theory and dielectric spectroscopic experiments for ω > 2.5 × 10 2 MHz 8 . As mentioned before the expected frequency for the collective oscillation of a biomolecule is around 0.1-1 THz, thus largely above the upper frequency threshold for important screening effects on SCieNtifiC REPORTS | (2018) 8:15748 | DOI:10.1038/s41598-018-33990-y electrodynamic fields. Collective phenomena are more probably expected in systems of resonant oscillators: for such a reason, a system of N identical biomolecules (oscillators) has been considered. Moreover, resonance of electric dipole oscillators, describing biomolecules, has been argued to be a necessary condition in order to activate long range dipole-dipole ( − R ij 3 ) electrodynamic interactions 5 . In our very simple model the force acting on each charge barycentre of the i-th electric dipole due to the j-th dipole is given by where E CED (r; R j ) is the value of the electric field in r generated by the j-th dipole whose center is in R j . According to the Classical Electrodynamics (CED), if we assume valid the dipole approximation, i.e. −  ‖ ‖ r r R j j , the expression for the electric field takes the form is the direction joining the center of dipole R j to r, p j (ω) is the Fourier Transform of the electric dipole moment of the j-th biomolecule in time domain and ε(ω) is the dielectric constant of the medium.
For the range of frequencies we explore (ω Ω ≈ 1 THz), the dielectric constant of an electrolytic aqueous solution can be assumed to be real ω ω  R I e ( ( )) m( ( ))   and approximatively constant ε WS (Ω) ≈ 3. Moreover both the intermolecular average distance R ij ≈ 10 3 Å and the characteristic linear dimensions r 0 ≈ 10 Å are much smaller than the characteristic wavelength of the electromagnetic field . This allows to assume that the electromagnetic field has the same value for both centers of charge of each biomolecule, i.e.
With these approximations the acceleration of the i-th dipole is directed along r i and due to the interaction with the j-th dipole reads as is the direction joining the electric dipoles and is a characteristic frequency describing the strength of the dipole-dipole interactions. Moreover is a geometrical factor depending of the orientation of the electric dipoles and r j (ω) is the Fourier Transform of r j (t).

Biological aqueous environment as thermal bath. This work is inspired by the request for observables
in real biological systems at molecular level that can detect the presence of long-range electrodynamics interactions among biomolecules. As all biomolecules in real biological environment are in aqueous solution, we have to take into account the presence of surrounding water molecules. Though recent studies reveal that the water in biological system can have a highly non trivial behaviour with respect to electrodynamic fields generated by the electric dipole of biomolecules [9][10][11][12][13] , in this article we will assume the surrounding water to play simply the role of a thermal bath. As a consequence of this, the presence of water molecules can be schematized via the introduction of a stochastic noise (thermal fluctuations) and a viscous friction term (dissipation) in the equation of motion for oscillating electric dipoles. In particular friction viscous forces are due to the aqueous surrounding medium considered as a homogeneous fluid with viscosity η w . We assume that the expression of the viscous force is given by Stokes' Law acting on each barycentre of electric charge (positive and negative) where  is the hydrodynamic radius of a typical biomolecule (~10Å).
SCieNtifiC REPORTS | (2018) 8:15748 | DOI:10.1038/s41598-018-33990-y From Eq. (8) it follows that the acceleration on the dipole length is given by On the other hand the stochastic forces are due to the collision of water molecules and freely moving ions on the biomolecules and they correspond to the realization of a thermal bath at temperature T. In particular these forces, acting directly on the charge barycentres of each biomolecules, can be described according to the following expression where ξ i (t) represents white noise whose characteristics along each Cartesian component α, β = x, y, z are given by The minus sign in the correlation term is due to the constrain we impose for thermal noise Such a condition does not take place in general for a real physical system but it has been implemented to provide a consistent realization of a stochastic system such that the center of mass of each molecule is fixed. With this prescription the stochastic force along the dipole direction is given by External forcing to produce out-of-thermal equilibrium conditions. In 5 it has been shown that long-range interactions among biomolecules can be present if the system of oscillating dipoles is maintained in out-of-thermal equilibrium. To achieve this goal a forcing term F NE,i (t) has been included in the equations of motion for the electric dipoles in order to ensure an external injection of energy. The explicit form of the force F NE,i (t) depends on the specific process that is chosen to inject energy into the system. In particular, a possible mechanism that has been used recently in THz spectroscopy experiments to detect collective giant oscillations in biomolecules, is the injection of energy in vibrational modes through the vibrational decay of the excited fluorochromes (i.e. fluorescent dye molecules) attached to each biomolecules 6 . This process can be represented choosing the following explicit form for the forcing term Equation of motion for the system of oscillating interacting dipoles. The equations of motion that describe the dynamics of the system with mutually oscillating dipoles are where all the biomolecules are assumed to be identical so that they all have the same characteristic frequencies ω i = ω 0 and Λ i = Λ.
In order to simplify the discussion we rescale the system according to that substituted in Eq. (18) give the following system of stochastic differential equations of first order where N is the total number of dipoles and 〈 〉  d is the average intermolecular distance in λ units. As a reference case in our simulations the parameters have been chosen to be = ∼ m 10, Z i = 1000, while the average intermolecular distance λ 〈 〉 = 〈 〉 = . × = . × −   d d m 1 6 10 Å 1 6 10 3 7 . The reason for choosing such a large value of Z is justified under the hypothesis that the surrounding water molecules participate to the effective dipole of each biomolecule and enhance it 5,9 . Therefore for the considered choice of parameters Ω . × − 2 3 10 ij 2 3 . Finally, in order to consider different cases with stronger interactions (corresponding to shorter average intermolecular distances, for instance) the coupling term is multiplied by a factor K > 0 with respect to the reference case just discussed, i.e.  This paper is intended as a first feasibility study for the detection of long-range electromagnetic interactions amongs biomolecules in aqueous solutions via a spectroscopic observable. In laboratory conditions the only parameter concerning molecule space configuration and orientation that one can easily control, is the biomolecules concentration, i.e. the intermolecular average distance. Therefore we investigate the emergence of collective behavior between biomolecules by varying the average distance among dipoles that corresponds in our model Eq.
(24) to vary K. On the other hand, in this work we do not investigate the role played by spatial correlation of the position and orientation of the dipoles in the appearance of a spectroscopic signature of the long-range electrodynamic dipole-dipole interactions. The study of dependence of spectroscopic observables on spatial correlations could be very interesting in this framework but we postpone to future work this investigation.
The parameter pul  can be estimated assuming that the energy injection on each biomolecule is due to the vibrational decay of a fluorescent dye. It is realistic 6   . Finally, the characteristic frequency for the energy transfer Ω pul is one of the most delicate parameters to be settled. This term accounts for the continuous injection of energy into the system, however the release must be done without perturbing too much the oscillating behavior, therefore we can assume that Ω Ω − 10 i pul 2 .

Numerical Results
The reported analyses have been done using a single system size (N = 50) and random initial conditions both for positions and velocities. However, similar results have been obtained for N = 100, 200 (not shown). The collective evolution of the population and in particular the level of coherence is usually characterized in terms of the macroscopic field where the modulus r 1 is an order parameter for the synchronization transition being one ( − N ( ) 1/2 ) for synchronous (asynchronous) states, while Φ is the phase of the macroscopic indicator 14,15 . However, in our case, the molecules are pivoted to the center of mass and cannot rotate: the effective degree of freedom of these objects consists in an elongation/shrinkage along the direction identified by the mutual distance between the two centers of charges. Therefore it is not possible to describe the movement of the dipole in terms of an oscillator rotating along the unit-circle via the identification of a time-dependent phase. The solution that we have adopted is to calculate the phase of the single molecule by using the inversion formulas In the following we present two set of parameters: a first one corresponding to the values discussed in Sec. II F and a second one, where we arbitrarily increase the thermal noise to investigate the robustness of the system and the emergent collective effects. The first set of parameters has been used to find the results shown in Figs (1-4). In this case the calculation of the order parameter r 1 does not lead to the identification of emergent (phase) synchronization in the system; in particular r 1 does not show any dependence on the coupling constant (see Fig. 1, panels (a) and (c)), as we would expect when the molecules are interacting with increasing strength. On the other hand if we calculate the order parameter usually employed to identify the emergence of 2- , we clearly observe a transition to cluster synchronization already for small coupling constants (see Fig. 1(c)). This is confirmed if we calculate the distribution of positions and velocities of the molecules (see Fig. 2, panels (a-i)). In particular, if we look at the phase space (x, v) it emerges clearly that the system splits naturally in 2 clusters and the distance between the clusters increases if the single elements are coupled. For stronger coupling more clusters emerge leading to a partial cluster synchronization. Moreover the probability distribution functions of the positions reveal an increase of the elongation towards values that turn out to be not realistic from a bio-physical point of view if the coupling constant is too big  K ( 1). Only the probability distribution profile of the velocities remains unchanged if the coupling constant and, therefore, the distance among the dipoles, is changed.
A better insight of the emergence of a collective behavior can be achieved if we consider the total dipole moment Due to the fact that the system is not deterministic and a white noise source is present into the differential equations, we have developed a method similar to the second-order Runge-Kutta one for solving numerically ordinary differential equations. In particular we have implemented the Heun method 16 in the Runge-Kutta algorithm as suggested in 17 , and we have used an integration time step 0.002 to perform the simulations. In addition to this, in order to compare the results for different coupling constants and for different strengths of the thermal ( ), it is possible to filter the low-frequency components of the spectrum just using the Fourier transform of the derivative. The low-frequency components that we want to filter out are related to the injected white noise that are not interesting for the scope of this work, i.e. finding a mark of emergent collective behavior.
Therefore we calculated the power spectrum of dP/dt to investigate the role played by the interactions among the dipoles to enhance a collective motion. While non-coupled dipoles show a peak at frequency ≈0.280 ± 0.006 (Fig. 3(a)), as soon as a small coupling is present in the system, the interactions among the dipoles get stabilized and a peak at lower frequency emerges already for K = 0.25 (Fig. 3(b)). However the eigenmodes emerging for small coupling are destroyed for bigger coupling, where the non-linearity in Eq. (24) prevents this self-organized behavior at small frequencies (Fig. 3(e)). The peak at small frequency emerges again for K ≥ 19 (Fig. 3(f)), but the intensity shown is smaller than before. On the other hand the main peak moves to higher frequencies for increasing coupling, but the intensity is more and more depressed (Fig. 3(g,h)).
A summary is presented in Fig. (4), where the peak intensity and the corresponding frequency values are given as a function of the coupling; while the secondary peak emerging at low frequency value for K ≥ 19 remains . We refer to the Supplementary Information (Sec. II) for the analysis of the synchronization level for this second set of parameters. While in absence of interactions (K = 0), the system shows a single pronounced peak at frequency ≈0.488 ± 0.006, once the interactions are active (K > 0), another peak arises at smaller frequency ≈0.263 ± 0.013. By increasing the value of K we observe an increase of the peak at lower frequency, to which corresponds a decreasing of the peak at higher frequency: a collective motion is enhanced due to interaction, while the motion corresponding to the non-connected situation is depressed (see Figs 5, panels (a-h) and 6(a)). On the other hand the position of the peak (i.e. the corresponding frequency value) does not change significantly if we increase the coupling constant (see Fig. 6(b)); the more evident increasing ratio for K > 20 is related to the fact that power spectra become richer and richer for higher coupling and secondary peaks arise. One of these secondary peaks (the main one) emerging at bigger coupling constant is also reported in Fig. 6 (panels (a,b)), and it is termed "Third Peak".
Finally, if we analyze in more detail the behavior of the first peak, related to the emergent collective motion, as a function of the coupling constant, it is possible to identify two different scales, once the figure is plotted in log-log scale (Fig. 6(c)). In particular, the different scales present for low coupling constant (K < 5) and for sufficiently strong coupling (K > 10) denote a transition between two different dynamical behaviors: the cross-over between two different regimes, from the one dominated by individual asynchronous behavior, to the one dominated by collective motion, with strongly interacting oscillators, is thus compatible with these two different scales.
If we now investigate the role of the thermal noise strength, we obtain a stochastic resonance effect 18 : the signal at low frequency (≈0.28 ± 0.09) can be boosted by adding white noise to the signal, which contains a wide spectrum of frequencies. The frequencies in the white noise spectrum corresponding to the original signal's frequencies resonate with each other, thus amplifying the original signal (i.e. the signal at low frequency) while not amplifying the rest of the white noise. Furthermore the signal-to-noise ratio is increased, while the added white noise is filtered out thanks to the band-pass filter that we have implemented calculating the power spectrum of dP/dt. In particular the low frequency peak, that corresponds in our case to the collective motion, is more visible for thermal noise strength Ψ ∼ = 0.03, to which corresponds a maximum in the peak high (see Fig. 7 panels (a,b)).
This peak is depressed for higher temperature and less likely to be revealed. On the other hand the peak at high frequency (≈0.56 ± 0.22), corresponding to the dynamics of isolated dipoles, can be also boosted by adding white noise into the system, but it does not decrease as significantly as the former one for higher temperatures, thus meaning that the single dipoles in this model are able to react to big level of noise, even though this is physically not plausible, since we would expect that dipoles will break up for high temperatures.

Discussion
Let us now comment about the physical meaning, and about the prospective relevance, of the results described in the previous Sections. As repeatedly stated, the present study was motivated by the need of finding an experimental strategy -complementary to the diffusion-based one already discussed in 3-5 -to detect the possible presence of electrodynamic attractive forces between biomolecules. Such a possibility emerges in the following framework. By pumping energy in the biomolecules of an aqueous solution, that is by keeping these molecules warmer than the solvent (out-of-thermal equilibrium), when the input energy rate exceeds a threshold value, then all, or almost all, the excess energy (that is, energy input minus energy losses due to dissipation) is channeled into the vibrational mode of the lowest frequency. In other words, the shape of the entire molecule is periodically deformed resulting in a "breathing" movement 6 . In doing so the biomolecules behave as microscopic antennas that absorb the electromagnetic radiation tuned at their "breathing" (mesoscopic) oscillation frequency. But antennas at the same time absorb and re-emit electromagnetic radiation, thus, according to a theoretical prediction, these antennas (biomolecules) can attractively interact at a large distance through their oscillating near-fields, and through the emitted electromagnetic radiation, provided that these oscillations are resonant and thus, take place at the same frequency 6 . The still open question is whether these electrodynamic interactions can be strong enough to be experimentally detectable, which would mean, in the positive case, to have some prospective biological relevance.
In the model of an aqueous solution of biomolecules that we have tackled, each individual molecule is assumed to be driven to an out-of-equilibrium mesoscopic vibrational mode which, in turn, excites an attractive electrodynamic force field associated with a −1/r 3 potential, where r is the intermolecular distance. By setting the parameters of the model to physically realistic values, we have numerically investigated the effect of varying the strength K of the mutual dipole-dipole electrodynamic interaction. The new observed phenomenon, shown in Section III, is the appearance of a collective behaviour involving all the molecules of the system, which is identified through K-dependent spectral features of a suitable macroscopic observable P(t), defined in Eq. (28). Furthermore, the analysis of the level of coherence in the system, done in terms of the standard Kuramoto order parameter r 1 and of the 2-cluster order parameter r 2 , shows that the collective dynamics cannot be simply traced back to synchronization. In particular, the order parameter r 1 measures the level of synchronization, while r 2 measures the degree of cluster synchrony in the system. Both order parameters oscillate irregularly thus implying collective irregular dynamics, irrespectively of the regular, periodic (breathing) behavior of the single dipole units (see also Fig. 4(a,b) in the Supplementary Information). Due to the combined effect of an external periodical forcing, of the presence of noise and of the interactions among biomolecules, we observed a self-organized formation of phase clusters characterized by different velocities and leading to (possibly) chaotic or quasiperiodic behavior at the macroscopic level. The emergence of nontrivial collective dynamics in systems composed of elements whose evolution is extremely simple, is a well-known phenomenon observed in complex systems. In particular chaotic irregular behavior emerging from regular unit dynamics has been seen in 19 , while the emergence of quasiperiodic motion has been shown in systems of oscillators 20 , neuron models 21 and rotators 22 .
From an experimental point of view, this means that by performing spectroscopic measurements at different concentrations of the solvated biomolecules we could detect the presence of electrodynamic intermolecular interactions. Varying the concentration C of the solution entails the variation of the average intermolecular distance 〈d〉 according to the relation 〈d〉 = C −1/3 . And varying C would be a practical way of experimentally changing the parameter K of the model. The variable P(t), represents the ensemble average of the projection of the dipole positions in the cartesian coordinates system: this is a spectroscopically measurable observable. Moreover, being related with the overall dipole moment of the solution, it can directly probe the emergence of a collective behaviour of the solvated molecules, collective behaviour which can only be driven by the presence of intermolecular interactions. A spectroscopic approach would thus entail a dichotomic, clear-cut, answer: if nothing would change in the absorption spectrum of the solutions at different concentrations, this would indicate that the solvated molecules do not interact at a (large) distance, to the contrary, concentration dependent spectral features would mean that the solvated molecules interact at a distance. In conclusion, the results reported in the present work outline a very promising experimental strategy -complementary to the diffusion-based one -to ascertain whether biomolecules can interact through long-range electrodynamic forces.