Signal metrics analysis of oscillatory patterns in bacterial multi-omic networks

Motivation: One of the branches of Systems Biology is focused on a deep understanding of underlying regulatory networks through the analysis of the biomolecules oscillations and their interplay. Synthetic Biology exploits gene or/and protein regulatory networks towards the design of oscillatory networks for producing useful compounds. Therefore, at different levels of application and for different purposes, the study of biomolecular oscillations can lead to different clues about the mechanisms underlying living cells. It is known that network-level interactions involve more than one type of biomolecule as well as biological processes operating at multiple omic levels. Combining network/pathway-level information with genetic information it is possible to describe well-understood or unknown bacterial mechanisms and organism-specific dynamics. Results: Network multi-omic integration has led to the discovery of interesting oscillatory signals. Following the methodologies used in signal processing and communication engineering, a new methodology is introduced to identify and quantify the extent of the multi-omic oscillations of the signal. New signal metrics are designed to allow further biotechnological explanations and provide important clues about the oscillatory nature of the pathways and their regulatory circuits. Our algorithms designed for the analysis of multi-omic signals are tested and validated on 11 different bacteria for thousands of multi-omic signals perturbed at the network level by different experimental conditions. Information on the order of genes, codon usage, gene expression, and protein molecular weight is integrated at three different functional levels. Oscillations show interesting evidence that network-level multi-omic signals present a synchronized response to perturbations and evolutionary relations along with taxa.


Introduction
An oscillating multi-omic network is a complex interlacing of interacting elements, which could be defined gene/protein oscillators. These elements are able to produce oscillations of a certain frequency potentially involving several and different cell processes [13]. It is proved that biomolecules, such as genes or proteins, could exhibit oscillatory behaviors. Moreover, if combined, they could generate oscillatory circuits observable on the temporal axis. The gene/protein oscillations are shown to be controlled by underlying regulatory networks and impact to different scales. [23,33]. In [19] a comprehensive review is reported, showing that, in general, living cells have a pervasive dynamic behaviour where the key transcription and regulatory factors oscillate on and off repeatedly even when cells are in steady states. These oscillations can be detected in circuits of genes/proteins involved in stress responses, signaling, and cell development [18]. Many biomolecules exhibit network-level metabolic interactions coordinated with cell growth, chromosome replication and cell division [39]. For example, the Min oscillation is discovered to be fundamental in the E.coli cell division [20]. In addition, biomolecule oscillations are found in very complex and multi-periodic signals. [1] generate a de-noised waveform from multiple significant frequencies, to provide oscillation statistics including signal metrics and multi-periodicity quantification. Furthemore, genetic circuits present some interesting intrinsic dynamics [14]; for example, in signalling pathways, they are responsive to feedback loops and show functional plasticity [18]. Recently, pseudotemporal estimations based on the level of mRNA and/or proteins have been introduced to detect oscillatory gene networks from single-snapshot experiments [7]. Unfortunately, even if the complex dynamics of bacterial processes could be predicted by keeping track of bacterial functional adaptation to single-snapshot perturbations (control vs treatment), the oscillatory dynamics have not yet been explored enough. Thus, a complete mapping of regulatory and control mechanisms is not yet known. Furthermore, the lack of experiments along temporal axis make difficult to recognise oscillating biomolecules and their circuital interactions. [26]. However, network-level synchronisation could be outlined by the hypotheses that every biomolecule in a network could interact with any other; this causes that the shared biomolecule oscillators synchronise the signals on common fluctuations. There are few examples in nature for which these assumptions are fully verified, such as circadian oscillations [12]. On the other hand, in synthetic biology, artificial oscillators, although often showing poor accuracy [25,15], are one of the most promising research fields. In particular, artificial oscillators allow the creation of genetic circuits focused on the execution of logical programming in living cells. The E.coli repressillator experiment represents a clear example of how genetic regulatory networks can be designed and implemented to perform new functions [3,9]. In our previous works [4,5], we investigated the E. coli response to ≈ 70 perturbations by monitoring network-level oscillation changes from controls to treatments and discovering that at the network level, there is another type of inter molecular multi-omic oscillation associated to each single pathway and experiment. This type of oscillation, as far as we know, has not yet been sufficiently investigated on fixed-time analysis. In particular, our multi-omic oscillations could be described as a multi-periodic signal given by the variation of interacting biomolecular multi-omics. This variation should be intended in terms of sequence low-hight alternations of multi-omic values [5]. To clarify these points, in Section 2.1 and Figure 1, the identification procedure of these multi-omic signal is described. With respect to our previous works, the number of organisms is extended to 11. Also, the experiment cardinality has increased to the order of thousands. Previously, we demonstrated how to measure the structural relations between the genomic and proteomic layers and how these led to oscillatory variations in response to perturbations. On the contrary, in this paper, taking advantage of signal theory and communication engineering, ad hoc metrics to better quantify network-level oscillatory features are designed and the algorithms for their computation are provided on an online repository. In detail, two new change point detection algorithms (CP D) [34,37] are introduced. These algorithms are capable of managing both the complexity of the variable amplitude of the multi-omic signal and the multiple periodicity. Further analyses on network-level synchronisations based on our novel signal metrics are provided. In particular, the analyses are focused on the interactions through different pathways of the same organism and modulated by different condition contrasts (CC ) (single-snapshot mRN A experiments) [22]. Through our approach, it is possible to recognise the oscillating networks, eventually evaluating if they are synchronised (periodically and simultaneously activated) or not synchronised, and how this feature changes across taxa. Integrated multi-omics are created from the following single omics: the codon usage [32], mRNA amount contrasts [22] and the protein molecular weight. Our final dataset is composed of 2.830.722 multi-omic signals from 11 different bacteria on thousands of environmental experiments. Network-level oscillatory variations are analysed with three functional levels of granularity from KEGG orthology [16]. The results confirmed and extended our previous findings, by showing that network-level multiomic oscillations exist in bacteria. Moreover, we found additional clues to support that the oscillatory networks are synchronised showing a combined dynamic response to perturbations. Furthermore, the comparisons between the various bacteria succeed in highlighting, in a completely innovative way, how network-level oscillations could reflect the effects of evolutionary pressure. Moreover, even if the maintenance of the gene order is not well understood [29], this research could give new clues to its meaning, underlining its dynamical relations with the proteomic layer under the evolutionary pressure [36].

Methods
In subsection 2.1 the procedure for multi-omic signal identification is described. In Figure 1 an overview of this task is provided. In detail, the mRNA condition contrasts and protein weight information are extracted, normalised and aligned with the codon usage information. In Figure 1 -Box (c), the multi-omic signal identification in condition con- Figure 1: The multi-omic information concerning one of the 11 bacteria is collected and data-integrated with respect to the j -th pathway and considering the k -th perturbation (indicated as the k -th experiment, or better condition contrast CC k ). As it is described in Section 2.1.1 and shown in Box (a) and Box (b), the multi-omic information is combined and grouped by molecular networks (KEGG pathways) and ordered with respect to the gene order information. The dashed lines between pathway nodes and genes indicate a 1:1 correspondence between the n selected pathway proteins and their related genes (∀i ∈ [1, n], p j,i ⇔ g j,i ). As shown in Box (c), the multi-omic values ( mv i j ) are represented as outputs of gene/protein oscillators at time t 0 ≤ k ≤ t f , which is the time when a single condition contrast is taken. Then, the ( mv i j ) are discretized into N quantisation levels (see also Section 2.1.2), and rotated of 90 • with respect to the multi-omic space. In this way, as it is shown in Box (d), a discrete multi-omic signal s i j is obtained, which could be eventually operon compressed (see Box (e) and also Section 2.1.3). The * 1 and * 2 asterisks represent the proximal and the distal positions of possible promoters or repressors neighbouring the operon g j,3 , g j, 4 .
trast CC at a fixed time k is indicated. From this single experiment, the collected multi-omic values mv are selected with respect to the KEGG pathway composition (Box (a)) and ordered considering the related DNA spatial positions (Box (b)). Then, they are figuratively rotated of 90 • and normalised generating a network-level multi-omic signal (Box (d-e)). As shown in Box (d), the multi-omics are discretised by considering N common levels of discretisation obtained through a between-organisms analysis (IOAC procedure) as defined in subsection 2.1.1.The identified signals have been found to be almost periodic, then in subsection 2.2, two novel algorithms are introduced for their analyses. Furthermore, by applying Algorithm 1 across the organisms and for all the multi-omic com-binations, the periodicity in the multi-omic signals is identified through an estimateθ of the search window where it can be found (Figure 2 -Box (a) and Section2.2.1). Next, another CPD algorithm (Algorithm 2) is applied to obtain two signal indices: oscs and osc k (Figure 2 -Box (b)). In particular, they are used to identify network-level oscillations for each pathway, for each experiment and for each multi-omic combination (Figure 4-Box (a)). Furthermore, the oscillations between networks and among the experiments, as described in the pipeline of Figure 4-Box (ab-c) are computed within and between organisms. Finally, this step is well described in Section 2.3.

Definition and normalisation
Given Pj,n as the j -th bacterial pathway of n proteins, the multi-omic signal sj is composed by a finite vector of n multi-omic values mvj associated to a subset of genes Gj,n = {gj,1, gj,2, . . . gj,n}. The gj,i are collected in Gj,n considering the exact correspondence 1:1 with respect to the proteins pj,i that compose the Pj,n = {pj,1, pj,2, . . . pj,n}. Each multi-omic value mvj,i is arranged on the multi-omic signal sj,i considering the relative position of its associated gene gj,i with respect to the origin of replication. In order to describe the signal we adopt the form sj [k] with the index k ∈ Z (such for example, the i -th multi-omic value of a signal is equal to sj[i] = mvj,i). The multi-omic values mvj,i are combined averaging their associated single-omics svj,i, ∀i ∈ Z . Since the svj,i are not defined in the same range, a normalisation is applied to make their values comparable. In particular the single-omics are three: (I) the mRNA condition contrasts (CC ) and (II) the molecular weights (MW) and (III) the codon adaptation index. MW and CC are normalised into the interval [0,1], that is the same range in which the codon adaptation index (CAI) is already defined [32]. In section 2.4 there is an accurate description of these sources.

Multi-omic discretisation through the IOAC procedure
In order to compare the signals between different organisms, an amplitude discretisation process was applied. The whole procedure is called: Inter-Organisms Amplitude Consensus (IOAC) and discretisation. The signal ) represents the c-th interval in which the mvj,i falls. However, the correct estimation of N for the discretisation follows a study on the single omic distributions trough the IOAC. In particular, the distribution of the svj,i on the whole genome and for each organism was investigated by means of the Anderson-Darling test (A-D Test, [28]). In the case of CAI, for the 91.67%, the A-D test rejects H0 with a significance of 0.05 . The MW and CC have the same significance with a percentage near to 100%. We conclude that the single omics do not follow a normal distribution. As a consequence, the optimal number of N bins is computed for nonnormal distributions applying the Doane's formula [38]. The N values are estimated between-organisms. In particular, N is fixed equal to 9 troughs a frequency based consensus (bold column in Supplementary Material Section 1). Consequently, in our set-up, the multi-omic signals generated for each experiment are discrete non-deterministic signals that represent gene-ordered multi-omic values that fall into 9 possible class intervals (from 0 to 8).

Operon compression
The extent of the multi-omic signal dataset is increased by their operon compressed versions. In this case, if a set of multi-omic values in a signal sj of length n are part of an operon in position r of length m, this set is defined as sj[r : r + m] = {mvj,r, mvj,r+1, mvj,r+m}, with m < n. In this case, each mvj,i follows its natural adjacent disposition on the DNA sequence. For this reason, we can represent the signal as a concatenation (indicated as ⊕) of the original signal with respect the operon: sj[1 : According to their natural regulatory functions, in order to apply the compression, the mvj,i that composes an operon could be seen as a single averaged value : sj[r] = fm(|mvi +mvi+1+. . .+mvm/m|). In this way we obtain additional signals with this shape: Obviously, the compression is applied more times if occurs, thus shortening the signal length.

Variable half-periods estimation with a change-point detector
In this section, we deal with multiple-periodicity signals characterised from different mRNA condition contrasts (see Section 2.4). These could represent more replications of the same experiment. Thus, due to the experimental intrinsic and extrinsic noise [35], it is very rare that these signals follow an ideal shape with a fixed periodicity; on the contrary, the periodicities are more variable making difficult the oscillation detection. After all, if there are oscillations, then their half-periods (from peak to lows or vice versa) are localised in windows of variable length. In order to obtain an estimation of half-periods, we introduce a novel localised and non-parametric change-point detection algorithm (CPD). The idea behind Algorithm 1 is based on the analysis of the median multi-omic variations estimating the window lengths θs in which the half-periods occur. A change-point is detected if and only if the median value of the previous variations ( m k ) is less or equal to the new coming multiomic variation (d k ) respecting the genes order. In general, the change point detectors look for changes in the statistical characteristics of the signal (i.e. the median values m k ), therefore considering the signal as a collection of different distributions arranged in adjacent windows [8]. In Figure 2 -Box (a) an atomic example of a half-period length estimation (θ) is shown. Algorithm 1 collects in a vector θ all the θs computed along the signal. Then, for each signal, the median values of θ are collected. Next, a commonθ is defined as the maximum between the median values of θ between all the organisms and for each multi-omic combination. In particular, theθs have a double functionality; they are indices of the different multi-omic interplay within and between organisms and they are halt condition parameters of Algorithm 2. In our case, we are focusing on the variable half-period lengths for all the pathways. It is discovered that all the organisms, on overall experiments and for each pathway, show a median θ of about 3.0 and an average comprised between 3.0 and 3.7 with a low standard deviation. Also, the max and min θs values are very similar between organisms. These statistics depend on the different multi-omic combinations (MOC) and on the presence of operon compression. The table of theθs is shown in the Supplementary Material Section 8, while in Section 2 the source code of Algorithm 1 is provided.
Algorithm 1 Multi-omic median window periodicity with CPD Require: A multi-omic signal: s[n] of length n and N the between organisms bins estimation In this section, we introduce algorithm 2, which is a multi-level change point detector, capable of detecting multi-omic variations in relation to the half-periods to which they belong. The algorithm returns in output oscs and osc k . In particular, oscs is an oscillation index which relates the length of the half-periods (conditioned by the multi-omic variations) with their signal amplitude. Instead, osc k is an index describing the relative length of the half-periods with respect to the signal length. The algorithm is divided into five steps, as they are summarised in Figure 2 -Box (b) and detailed as follows: I) Collect multi-omic adjacent change-points: In this step algorithm 2 collects the adjacent multi-omics change-points on the signal s tracking the adjacent multi-omic class variations am as described in Equation 1:

II) Collect adjacent multi-omic variations:
For each signal s, algorithm 2 collects the adjacent multi-omic variations amv as the absolute difference between two adjacent multi-omic values, as described in Equation 2: .
The calculation of mvq and mvl is the central point for the development of our change point detector, therefore Step III has been described in detail in the pseudocode of Algorithm 2 and is illustrated in the example of Figure 2 for an halt condition. IV) Computation of oscs: The two vectors mvq and the mvl have the same length d, where at most d ≤ |s| − 1. The first one traces the quantified multi-omic variation for each half-period. The second one traces the window lengths in which the variations are computed. Thus, the oscillation index is defined in Equation  5 as the product of the two vectors with respect to the sum of mvl for the N bins: The : For each signal, this index represents a relation between the length of s and the number of periods along the signal, described as the cardinality of mvl. As we will see, osc k ∈ [0, 1] remains defined in a certain interval and describes some interesting relations in the analysis of the pathways phase synchronisations (Section 3). In the Section 5 of Supplementary Material we provide the source code related to Algorithm 2 in order to compute oscs and osc k .

Robustness and sensitivity analysis of Algorithm 2
The robustness of Algorithm 2 was tested by defining two types of perturbations. Without loss of generality, we assumed that the perturbations are defined by a random distribution with zero mean and unit variance. As a consequence, the first type of perturbation applied to the discrete signal s is a stochastic additive noise. We decided to add the 5% of the generated noise, in the following way: s[i] + (N (0, 1) * 0.5) ∀i ∈ [1, n]. The second type of perturbation consists of random shuffling the elements of the original signal, thus testing the importance of the information deriving from the gene order. The t-test p-value of the obtained oscs on random shuffled distributions is equal to 0.004385, while on random additive noise distributions it is less than 2.2e − 16. We selected a subset of original signals with at least 70% of significant oscillating multi-omics. This means that we selected only the signals with oscillation index oscs ≥ φ, with φ = 0.7. Then, we perturbed this subset of signals and we newly computed the oscs. In Figure 3, the PDF of the original signal oscillation index oscs (solid line) against the perturbed ones (dashed lines) are shown. In particular, we can observe that the oscs of the original signals remains into the interval from 0.7 to 1, while the noisy and shuffled signals intercept an interval from 0.4 to 1, thickening the area of interest ( Figure 3 - * 2 , * 3 ) Algorithm 2 : Oscillation indices: osc s and osc k . The procedure is divided into 5 steps as described in Section 2.2.2. Here, in pseudocode, Step III in relation with the other steps is shown.
Require: A multi-omic signal: s[n] of length n and estimatedθ function compute-mvq-and-mvl (s[n],θ) am ← Step I ⊲ Collect multi-omic adjacent change-points. amv ← Step II ⊲ Collect adjacent multi-omic variations. to lower values than those defining the original area ( Figure 3 - * 1 ). Note that Algorithm 1 and Algorithm 2 have linear complexity O(n) over the signal length n. As expected, small variations in multi-omic values or a random arrangement clearly lower the oscillation index oscs.  The multi-omic signals were integrated for every organism (∀O) and for every Pj following the combinations of single omic values svj,i with and without operon compression. The interplay between these omic layers is described by [2]. In particular, these multi-omic values ( mvj,i) are combinations of dynamic svj,i and static svj,i or only a combination of static svj,i. The static svj,i are the CAI (genomic layer) and the MW (proteomic layer).

Multi-omic sources
The svj,i of the CAI was computed as described by [32]. The MW was computed considering the molecular weight of the amino-acidic composition of each protein in Pj . These values are called static because they do not change when perturbations occur. The dynamic svj,i were represented as 'condition contrast' (CC ) and represented the mRNA expression changes between microarray experiments. From a certain point of view, the CC is a glue between the two static layers. The data for the condition contrast were downloaded from the COLOMBOS v3.0 dataset and they were already normalised within and between experiments by the creators of the dataset [22]. For this reason, through the CC was possible to compare the between organisms signals s∀O based on multi-platform experimental setups without losing generalisation. The four multi-omic combinations considered in this work were: The molecular networks pathway information and their KEGG Orthology are extracted from KEGG through a REST service [16,21,17]. The information about the operons localisation was determined through the OperonDB dataset [24]. The order of the genes was obtained through the NCBI dataset [6] and aligned to the KEGG microbial genome information. The whole dataset, with their respective labels, is provided in an R data format in the Supplementary Materials Section 6. The phylogenetic tree is reconstructed considering the NCBI taxonomy dataset [10].

Results and Discussion
I) Multi-omic signal oscillations: The mean absolute error (MAE) is computed in order to quantify the distance of oscs obtained from the original signals and the noised/perturbed ones (see also Supplementary Material -Section 4) [31] . The two types of perturbations adopted here are the same described in paragraph 2.2.3. MAE is computed for all the s ∈ O ( see Section 2.4) with oscs ≥ φ with φ = 0.7. Algorithm 2 computes the oscs of the original signal and of the noisy one in the same way. A cutoff is applied to the wavelength, considering that the biological meaning of sequences makes sense for signals with n ≥ 6. Nevertheless, very long signals present low oscillation indices with respect to median wavelengths, thus, respecting the short memory property of [30]. With these constraints, we are still considering the 95% of the pathways and 63% of pathways with operon compression. The oscillation index oscs is computed on signals with and without operon compression considering all the possible MOC separately. The MAE between the original signals and the perturbed ones is shown in Table 1. As it is shown, the distance between the noisy signals and the original ones is more pronounced considering the combination of CAI, molecular weight and condition contrasts (bold cells). Instead, the distance between the random shuffled signals and the original ones is more pronounced in the combination of CAI and MW (bold cells). In this case, the involved omics are both static and deeply related to the gene order. In these analyses, it is highlighted that the multi-omic signals preserve the oscillation behaviours, described by our oscillation indices, representing effectively recurrent patterns present in nature. As it is shown also in Figure 3), slight or massive multi-omic variations lower dramatically the oscillation index oscs. II) Network-level synchronisations: The synchronised signals, for each experiment, could give a meaningful picture of the network-level interplay of multi-omics. With the methodology described in section 2.3, important clues to support the hypothesis that there are groups of pathways/oscillatory networks in synchronisation are provided. For example, in Figure 4, the most frequent network-level interplay based on multi-omic signals CAI-MW-CC for the functional feature of KO Level 2 is shown. The intersections reported in Figure 4 involve the pathway signals of all 11 bacteria, shaping the general behavior of the dynamics of oscillations at the network level and outlining a similar response to several experimental perturbations. According to [11], if we look at, for example, Figure 5 box The distributions of the osc k were studied by separating the synchronised signals (gray boxplots) to those with a lower oscillation behaviour (white boxplots). In this case, the oscillation index osc k of the synchronised signals is lower than in those not synchronised even if they show a higher oscillation index oscs. It is evident that the two distributions are separable and preserved along with the organisms, except in rare cases where some few oscillatory networks seem to overlap the not synchronised ones.
In these particular cases, we observed a higher osc k for these functional classes:   Table 1 the mean absolute error (MAE) between the original signal oscillation indices (osc s ≥ φ, with φ = 0.7) and the perturbed ones are shown. The additive noise and the random shuffle are the two perturbations listed as noise and shuffle. The distances are computed for each MOC. In Signals setsize, the size of the signals with the in-text described constraints is reported. Despite the set-size that remains equal for each organism, every pathway and every MOC, the number of multi-omic operon compressed signals is different from the original ones because some pathways during the operon compression became too short to be considered as biological sequences.
dence could allow us to assume that network level synchronisations have a significant oscillation in amplitude with longer periods due to the underlying synchronisation mechanism that can be derived from the regulatory and control circuits. Conversely, a high osc k could indicate that networks show a rapid biological response to external stress. However, these hypotheses will have to be explored in future research work.

Conclusion
In this article, two new signal metrics have been introduced to study multi-omic oscillations at the network level and on single-snapshot experiments, defined in particular as condition contrasts. From the analysis of these metrics, it is possible to provide interesting clues about the characteristics of the signal proving that multi-omic network-level oscillations exist in nature. Furthermore, clear evidence has been provided that these oscillations could show a synchronised interaction in response to perturbations. Multi-omic signal analyses have been extended to multiple organisms and related to their phylogenetic tree to provide better comparisons. Algorithmic methodologies are provided and accompanied by a tool on supplementary material and on the online repository. This work could be useful in the fields of synthetic biology and systems biology with the goal of mapping the organism regulation and control circuits, for example, in case of lack of time series experiments. Furthermore, there is a growing amount of whole-genome and longitudinal data and these metrics answer to the need to detect complex patterns of changes.  It is plotted the PDF of the original multi-omic signals with osc s ≥ φ with φ=0.7 (solid line). As it is proved, the underlying area * 1 is comprised into the interval of osc s from 0.7 to the upper bound of 1.0. Then, this subset of original signals (osc s ≥ 0.7) is perturbed in two ways. The PDF of the original signals perturbed with noise are shown as dashed lines and the underlying area is indicated with * 2 . Those random shuffled are shown in 3 dots dashed lines and the underlying area is indicated with * 3 . It is possible to observe that, when the original signals are perturbed, the PDFs area * 3 and * 2 move mostly on osc s values comprised between 0.4 and 0.7. Thus, the perturbed signals lower the osc s proving the Algorithm 2 correctness, in terms of robustness and sensitivity analysis. The osc s is computed for all the signals with a length of at least 6.