teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The history of mankind cannot be written without mentioning the impact that viruses have had on it [START_REF] Diamond | Guns, Germs, and Steel[END_REF]181]. As put forward by Flint et al. [90], viruses has inuenced on human history as much as any war, natural disaster or invention. The inuence of viruses on human lives is not recent [221]. It is nowadays known that as far back as AD 600, humans carried the variola virus that produced smallpox [170], which was possibly circulating in humans since around the fall of the Western Roman Empire, when many peoples were migrating across Eurasia. It has been recently shown that measless virus Institute for Cross-Disciplinary Physics and Complex Systems (IFISC, UIB-CSIC), Campus Universitat de les Illes Balears E-07122, Palma de Mallorca, Spain. E-mail: estrada@isc.uibcsic.es might have jumped to people in the rst millennium BC [START_REF] Düx | Measles virus and rinderpest virus divergence dated to the sixth century BCE[END_REF], and hepatitis B had been infecting humans since the Bronze Age, 5,000 years ago [169].

In the XX century along, smallpox killed over 300 million people worldwide, which is more than twice the estimated number of deaths from all wars in that centuryat least 108 million people [181]. In 1918In -1919 the so-called Spanish inuenza pandemic caused approximately 50 million deaths worldwide [168,229]. More recently, pandemics like the one produced by HIV-1 has infected more than 60 million people producing more than 25 million casualties [213]. Viruses have been producing pandemics regularly through history. In the review paper Pandemics through history Piret and Boivin [186] compiled several of the most relevant pandemics that mankind has suered through history. We have picked those which are produced by viruses and represented them graphically in Fig. 1. Since the end of the XIX century to the end of the XX one there were ve major pandemics spaced at 39, 29, 11 and 13 years, respectively. However, in the rst 20 years of the XXI century humans have suered four major pandemics which are spaced at 7, 6 and 4 years, respectively. The four pandemics of the XXI century described in Fig. 1 are of zoonotic origin.

Zoonoses are infections transmitted between species, in particular from other animals to humans [136]. The question about whether the frequency of zoonotic infections in recent years is increasing or not is out the scope of this chapter and need a detailed analysis of a wider literature. However, by considering a few of the parameters that govern the capacity of a virus to infect a human population we can learn a few hints about what is happening and will happen in the future. The dramatic increase in traveling in recent years has supposed a signicant increment in the number of contacts between dierent parts of the world population [START_REF] Brockmann | The hidden geometry of complex, networkdriven contagion phenomena[END_REF]. Many infectious diseases are transmitted through vectors, as it was rst predicted by Cuban physician Carlos J. Finlay [START_REF] Finlay | The mosquito hypothetically considered as an agent in the transmission of yellow fever poison[END_REF] and corroborated by the Reed Commission for the case of yellow fever transmitted by mosquitoes [START_REF] Chaves-Carballo | Carlos Finlay and yellow fever: triumph over adversity[END_REF]. Increase of global travels and transport of goods across the world, may produce that vectors harboring a given virus can escape their natural habitat, conquering new ones, where the virus can spread [START_REF] Flint | pathogenesis and control[END_REF]. In the same direction, the changes in climate, e.g., temperature and humidity, can make that the natural habitat of a given vector, e.g., a mosquito, which before was constrained to certain tropical or subtropical regions, now expands its frontiers to wider world regions. Another factor facilitating the virus transmission is the increase in population density [175]. With a greater concentration of humans in cities, the increase in population density is an ever growing factor that facilitates virus spread. Finally, we would like to mention the increasing damages of ecosystems due to human activity [127]. According to Quammen [192] we are interacting with wild animals and disrupting the ecosystems they inhabit all to an unprecedented degree. This ecosystems disruption may imply that we are more exposed to zoonotic infections as it is evidenced by the large number of such infections since the middle of the XX century, e.g., Macupo virus (1961), Marburg fever (1967), Lassa fever (1969), Ebola (1976), HIV-1 (1981), Sin Nombre virus (1993), Hendra virus (1994), avian u (1997), Nipah virus (1998), West Nile virus (1999) SARS (2003), MERS (2015), SARS-CoV-2 (2019).

Viruses also aect other animals, particularly those of economic value. As an example the outbreak of foot-and-mouth disease in the U.K. in 2001 resulted in the sacricing of over two million heads of livestock and economic losses estimated at ¿6. 5 billion [212]. Plants of human value are also the target of viruses [120,180,222].

Aphids (aphidae) [START_REF] Blackman | Aphids on the world's crops: an identication and information guide[END_REF]238] alone are reported to transmit about 50% of insect-borne plant viruses (approximately 275 virus species). From the economic point of view this virus transmission by aphid represents global losses estimated on tens of millions to billions US$ of yield loss per annum [START_REF] Blackman | Aphids on the world's crops: an identication and information guide[END_REF][START_REF] Dedryver | The conicting relationships between aphids and men: a review of aphid damage and control strategies[END_REF]129]. In the UK alone the damage on cereals made by aphids has been estimated to be around 60-120 million pounds annually [228]. Due to climate change pests and pathogens have shifted poleward by 2.7 ± 0.8 km/yr since 1960. Deutsch et al. [START_REF] Deutsch | Increase in crop losses to insect pests in a warming climate[END_REF] estimated that global yield losses of rice, maize and wheat grains are projected to increase in the range of 10 to 25% per degree of global mean surface warming. Thus, in a projected scenario of 2°C-warmer climate the mean increase in yield losses owing only to pest pressure extend to 59, 92, and 62 metric megatons per year for wheat, rice and maize, respectively [START_REF] Deutsch | Increase in crop losses to insect pests in a warming climate[END_REF].

What are our weapons to ght existing and new viruses to come? At individual level humans have their own genetic and immune systems which are the rst biological contention barrier against viruses [173]. However, our systems are not necessarily prepared for ghting against quickly evolving viruses [START_REF] Ahmed | Protective immunity and susceptibility to infectious diseases: lessons from the 1918 inuenza pandemic[END_REF]. Decreasing nongenetic risk factors is an important way to be prepared, as having a more healthy system always improve our defenses [START_REF] Biondi | Psychological stress, neuroimmunomodulation, and susceptibility to infectious diseases in animals and man: a review[END_REF][START_REF] Flint | pathogenesis and control[END_REF]. But without any doubt, the most important of all mankind weapons against future pandemics is science. Virologic and epidemiologic surveillance [111] may be intensied to detect the potentially emerging new viruses and to establish procedures to shorten the periods between the beginning of an epidemic and its detection. New strategies to create vaccines in short time, as it has been successfully done for the case of SARS-CoV-2 [START_REF] Golob | SARS-CoV-2 vaccines: a triumph of science and collaboration[END_REF]223,240], is also of vital importance. Obviously, there are many scientic disciplines involved in this scientic endeavor: virology, molecular biology, anatomy, physiology, epidemiology, and more. On the grounds, supporting all of them, is the emerging science of complex systems [START_REF] Boccara | Modeling complex systems[END_REF], which has been claimed to be important for public health studies in general [203]. This integrative science deals with the study of those systems formed by a large collection of interacting parts which usually show self-organization with non-trivial structures and dynamics. As we will show in this chapter an overarching concept for the study of virus-host interactions from a holistic view is that of networks [START_REF] Estrada | The structure of complex networks: theory and applications[END_REF]141,176]. Therefore, we will cover here network-theoretic studies ranging from the molecular biology of viruses, cellularlevel virus-host interactions, organism-level eects of virus infections, and society-scale transmission and impact of viruses. Due to the importance of the current pandemic produced by SARS-CoV-2, which is known as COVID-19, we will focus our studies on this virus/epidemic. There are some topics, which although important, will not be covered in this Chapter. One of then is the use of network theory for studying RNA structure. The second is about network vaccinology. The rst topic was reviewed by Kim et al. [131] in a book chapter in 2013 and the second topic was previously partially covered in a review by the current author [START_REF] Estrada | COVID-19 and SARS-CoV-2. modeling the present, looking at the future[END_REF], particularly for epitope vaccines.

Viruses

Viruses are replicating organisms which are the most numerous of all biological entities in the world [START_REF] Breitbart | Here a virus, there a virus, everywhere the same virus?[END_REF][START_REF] Brüssow | The not so universal tree of life or the place of viruses in the living world[END_REF]. They are structurally very simple, appearing in various shapes, such as helical, icosahedral or in more complex ones (see Fig. 2) [200]. Because viruses lack a translational apparatus, they must replicate only inside other living organisms.

The variety of organisms infected by viruses is very wide and include Eukarya, Archaea, Bacteria [START_REF] Edwards | Viral metagenomics[END_REF]142] and even other viruses [139]. Most viruses which are known have diameters between 20 and 300 nm (1nm is 10 -9 m) [START_REF] Boºi£ | Statistical analysis of sizes and shapes of virus capsids and their resulting elastic properties[END_REF]. Only giant viruses like mimivirus are bigger reaching a capsid diameter of 400 nm and having protein laments measuring 100 nm at its surface [START_REF] Abergel | Giant viruses[END_REF][START_REF] Brandes | Giant virusesbig surprises[END_REF]. Some loviruses, which have lament-like shapes, can have length of up to 1400 nm, but with diameters of only about 80 nm [START_REF] Feldmann | Classication, structure, and replication of loviruses[END_REF]. All viruses have their genome encapsulated within a capsid [266]. It is a shell formed by several protein subunits known as capsomers (see Fig. 3). The capsid protects the viral genome from outside. This protective protein coat could be exposed to the exterior in naked viral forms, or can be covered by an outer lipid bilayer known as a viral envelope [253]. This membrane includes a number of membrane proteins and it is located above another proteinaceous coat, which is known as the matrix. Some viruses may have glycoproteins known as spikes either directly over the capsomers or on the envelop of the virus [146,200].

From a genomic point of view virus can have single-or double-stranded DNA or single-or double-stranded RNA [152]. They can be transcribed via positive sense, negative sense, or ambisense transcription mechanisms. A classication of viruses in seven Fig. 3: Main structural characteristics of naked and enveloped viruses.

classes according to the nature of their genetic material was provided by Baltimore [START_REF] Baltimore | Expression of animal virus genomes[END_REF] (see also [134]) and it is resumed in Fig. 4 SARS CoV-2 is a Betacoronavirus, an enveloped, single-stranded RNA virus of positive (+) sense (Class IV of Baltimore classication), it is of a spherical shape with a diameter of approximately 100 nm, a volume of approximately 10 6 nm 3 , and genome size of about 30 kilo bases, corresponding to the largest RNA genomes among viruses [START_REF] Bar-On | Science forum: SARS-CoV-2 (COVID-19) by the numbers[END_REF]201,239].

The replication process of SARS CoV-2 is illustrated schematically in Fig. 5, which will be used here as a guide for explaining the process [START_REF] Astuti | Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2): An overview of viral structure and host response[END_REF][START_REF] Glebov | Understanding SARS-CoV-2 endocytosis for COVID-19 drug repurposing[END_REF]215,243]. The replication starts when SARS CoV-2 enters human cells by means of receptor-mediated endocytosis mechanism, in which its spike glycoprotein S, previously preprocessed by the transmembrane serine protease 2 (TMPRSS2), binds to the host receptor Angiotensin-Converting Enzyme-2 (ACE2) [START_REF] Abadias | Fractional-order susceptible-infected model: denition and applications to the study of COVID-19 main protease[END_REF]. In the interior of the human cell, the virus enters in the form of an endosome [START_REF] Abergel | Giant viruses[END_REF], where cathepsin L, an endosomal acid protease, activates the spike protein by cleaving the protein into S1 and S2. At this point, the S2 protein fuses the membrane of the virus with that of the endosome [START_REF] Ahammed | Estimation of novel coronavirus (COVID-19) reproduction number and case fatality rate: A systematic review and meta-analysis[END_REF], which results in the release of the viral material, RNA and proteins, into the human cell cytoplasm. Now, viral mRNA is translated (4) by using the human cell ribosome to form two viral replicase polyproteins: pp1a and pp1ab [START_REF] Akbar | Extracellular vesicles in metabolic disease[END_REF]. These two polyproteins are then cleaved by the Papain-like protease and 3C-like (main) protease, resulting in 16 non-structural proteins (Nsp) [START_REF] Alene | Serial interval and incubation period of COVID-19: a systematic review and meta-analysis[END_REF]. A replicase-transcriptase complex is then assembled with some of these proteins and the viral RNA [START_REF] Aleta | Modelling the impact of testing, contact tracing and household quarantine on second waves of COVID-19[END_REF]. It drives the production of subgenomic RNA(+) through transcription [START_REF] Moreno | Evaluation of the potential incidence of COVID-19 and eectiveness of containment measures in Spain: a data-driven approach[END_REF] and of pre-genomic RNAss(-) (8a) which is replicated to genomic RNAss(+) (8b). The subgenomic RNA(+) is now translated into structural proteins [START_REF] Alibhai | Cellular senescence contributes to age-dependent changes in circulating extracellular vesicle cargo and function[END_REF], such as the spike S protein, envelope proteins (E), membrane proteins (M), nucleocapsid proteins (N), and several Open reading frame proteins (Orf ). All of them enter the endoplasmic reticulum [START_REF] Alon | Finding and counting given length cycles[END_REF] where the nucleoprotein complex is formed between the nucleocapsid protein and the (+) strand genomic RNA. Finally, the assembly of all proteins and RNA into a new virus particle is carried out in the Golgi apparatus of the human cell [START_REF] Althaus | Estimating the reproduction number of ebola virus (EBOV) during the 2014 outbreak in west africa[END_REF]. The virus is then expelled from the cell via exocytosis [START_REF] Anand | Coronavirus main proteinase (3CLpro) structure: basis for design of anti-SARS drugs[END_REF] and starts its maturation to start a cycle again [START_REF] Anderson | Infectious diseases of humans: dynamics and control[END_REF].

Virus infections in the context of complex systems

The study of viruses is usually presented as a set of concentric circles [START_REF] Flint | pathogenesis and control[END_REF]. A rst circle may be provided by the detailed analysis of the genome and of the structures of viral proteins through molecular biology techniques [152]. It plays a fundamental role to understanding the biochemical mechanisms of virus-host interactions, which also overlaps with the system biology view of the problem [124,135]. A second circle, which partially includes the previous one, is provided by the understanding of viral replication mechanism and the analysis of how its infection aects host cells [173]. Increasing consecutive circles explain viral pathogenesis as infected cells are perturbed through virus-host protein-protein interactions, altering cell biochemistry and aecting the tissues in which infected cells are embedded into. Viral pathogenesis extends its circles through these tissues to larger anatomical systems, i.e., organs, which will malfunction as a consequence of the viral infection, then disturbing the physiology of the host as a whole [198]. A fundamental part of virus survivability consists in the transmission of viral populations from an infected individual to susceptible, not infected, ones. The analysis of an epidemics outbreak via inter-personal contagion mechanisms is the subject of study of epidemiology [START_REF] Giesecke | Modern infectious disease epidemiology[END_REF]125], which completes the main overarching circles of specialties (see Fig. 6) needed to understand viral infections in a host population. These specialtiesmolecular biology, virology, anatomy, physiology, epidemiology typically study separate parts of the complex mechanisms of viral infection of a host population. They need to cover huge spatial and temporal scales ranging from about 1nm ( 10-9 m) to several thousand kilometers (∼ 10 6 m), and from femptoseconds (10 -12 s) to months/years (∼ 10 7 s). Fig. 6: The size scales overarching the complex system of viruses infecting a host population.

We can then consider the global system of virus infecting a host population as a complex system . It is composed by several parts which are also complex systems in their own. Some characteristic features of complex systems are clearly observed here [START_REF] Boccara | Modeling complex systems[END_REF]. For instance, It a nested system which is formed by several components that themselves are complex systems; the system as a whole and its component parts are formed by many local interactions giving rise to networks of interactions between dierent elements of the system; the system as a whole and many of its components have memory, such that previous states in their dynamical evolution may have an inuence on future states;

the system is open, exhibits critical transitions, emergent properties, and contains feedback loops, such that the eects of the component of the system are fed back to that component altering its own state.

The nested nature of the system and its networked structure are graphically illustrated in Fig. 7. They can be considered as the main responsible of many of the dynamical and functional properties of the system as a whole. For instance, the patterns of infection of a given virus in a host population can be of dierent nature, producing either acute infections characterized by short reproductive cycles which produce many virus progeny, or persistent or latent infections, where the infected host survives for extended times such that faster viral reproduction does not provide selective advantage [173].

These patterns are conditioned by the own viral population dynamics. That is, let N (t) be the population of the virus at a given time. Then, the time evolution of this population can be described mathematically by the logistic equation (see p. 139 in [START_REF] Flint | pathogenesis and control[END_REF]):

dN (t) dt = rN (t) [1 -N (t) /K] , (1)
where r is the dierence between the average birth and death rates of the virus, respectively, and K is the carrying capacity of the environment (e.g., a single cell, an individual, the whole host population) which denes the upper bound of the population size of a virus that can be sustained by that specic environment. Then, the pattern of infection is very much determined by the magnitude of r and K. When r is relatively high and K is relatively low, the virus reproduction emphasizes high growth rate, which are characteristics of acute infections with extensive viral spread, such as in inuenza virus, norovirus, and rhinovirus [START_REF] Flint | pathogenesis and control[END_REF]. When r is relatively low and K is relatively high, the pattern of infection often correspond to those of persistent or latent infections like in Epstein-Barr virus infections [START_REF] Flint | pathogenesis and control[END_REF].

The pattern of infection also determines very much the evolution of infected hosts.

For instance, an individual who has suered an acute infection will be immune to challenge a reinfection by the same virus. However, some acute infections occur repeatedly on the same host due to the high capacity of structural proteins of these viruses to tolerate amino acid substitutions, i.e., structural plasticity, and/or their capacity to change surface proteins in response to antibody selection, i.e., antigenic variation [START_REF] Flint | pathogenesis and control[END_REF]173,198]. Here we observe feedback loops at dierent scales of the complex system of virus infecting a host where virus reproduction determines the infection pattern, which generates some kind of host immune response. Such immune response may trigger changes on surface proteins of the virus, which then may develop reinfections of the same host. Fig. 7: The networks covering the dierent scales of the complex system of viruses infecting a host population.

Networks

In this section we dene the most important concepts about networks which are then used in this chapter. Most of the concepts can be found in [START_REF] Estrada | The structure of complex networks: theory and applications[END_REF] where more references can also be found.

Denition 1 A network is a (weighted) graph G = (V, E, W, f ) where V is the set of vertices or nodes, which represent the entities of a complex system, E ⊆ V ⊗ V is the set of edges where for v, w ∈ V , (vw) ∈ E if there is a kind of relation or interaction between the two nodes, W is a set of weights assigned to the edges by the surjective mapping f : E → W . If w (e) = 1 for all edges e ∈ E the graph is said to be unweighted, otherwise it is a weighted graph. If for all (vw) ∈ E implies that also (wv) ∈ E, the graph is undirected, otherwise it is directed. Typically we will represent the graphs as G = (V, E), which would include the weights on the edges set.

Remark 1 Two nodes v and w are said to be adjacent in an undirected graph G if and only if (vw) ∈ E. In the case of directed graphs there are several denitions of adjacency depending on the author. Here we will prefer to use the term v points to w (equivalently, w is pointed from v) if there is an arrow from v to w.

Example 1 In Fig. 8 we illustrate four networks representing: (a) the spike protein of SARS-CoV-2, (b) the interactions between SARS-CoV-2 proteins with human proteins, (c) a social contact network of inhabitants of Haslemere, and (d) the airport passengers transport network of USA. In the case of the SARS-CoV-2 spike protein (a), the nodes represent the α-carbon atoms of each residue in the protein and two nodes are connected if the distance between such α-carbons is less than 7Å as determined by

x-ray crystallography. The crystallographic structure is deposited in the protein databank (PDB) with code 6VXX. The colors corresponds to the three dierent chains of this protein. In the protein-protein interaction network (b) of SARS-CoV-2 and human proteins, each node corresponds to a protein, red dots are for SARS-CoV-2 proteins and blue ones for human proteins. A SARS-CoV-2 protein is connected to those human proteins to which it is found experimentally to interact with. Some other human proteins, which are at two steps of separation from the SARS-CoV-2 ones, are also given.

The social network in (c) corresponds to 468 individuals (gray nodes) with 1,257 social links (blue edges) weighted by 1,616 daily contacts (edge thickness). A node marked in red was used by Firth et al. [89] as a single starting infector for simulations of a disease propagation on that network. Finally, in (d) the nodes correspond to airports in the USA and two nodes are connected if there is at least one ight between the two airports. In this case the edges of the network can be weighted by the number of passengers that ight from one airport to another in a given period of time.

Denition 2 The degree of the node v ∈ V , kv, in an undirected graph is the number of nodes adjacent to v. In a directed graph we call in-degree of v, k in v the number of nodes pointing to v, and its out-degree k out v the number of nodes to which v points out.

Denition 3 Two graphs G and H are isomorphic if there is a one-to-one correspondence between the nodes of G and those of H, such as the number of edges joining each pair of nodes in G is equal to that joining the corresponding pair of nodes in H. If the graphs are directed the edges must coincide not only in number but also in direction.

If the graph are weighted the edges must coincide also in their weights. Example 2 The two graphs two graphs in Fig. 9 are isomorphic as there is a one-toone correspondence between the nodes of the graph in a) and those of b), such as the number of edges joining each pair of nodes in the rst is equal to that joining the corresponding pair of nodes in the second. They correspond to two dierent drawings of the same protein residue network. Denition 4 The graph S = V ′ , E ′ is a subgraph of a graph G = (V, E) if and only if V ′ ⊆ V and E ′ ⊆ E. A particular kind of subgraph is the clique, which is a maximal complete subgraph of a graph. A complete graph is the one in which every pair of nodes are connected. Denition 5 A (directed)

walk of length L from v 1 to v L+1 is any sequence of (not necessarily dierent) nodes v 1 , • • • v L , v L+1 such that for each i = 1, 2, • • • , L + 1 there is link from v i to v i+1 . A walk is closed (CW) if v 1 = v L+1 . A particular kind of walk
is the path of length L, which is a walk of length L in which all the nodes (and all the A cycle is a closed walk in which all the edges and all the nodes (except the rst and last) are distinct. The girth of the graph is the size (number of nodes) of the minimum cycle in the graph. A cycle is chordless if there are no edges connecting pairs of vertices of the cycle except for those edges that itself belong to the cycle. Denition 6 In an undirected graph the shortest path distance d (v, w) = dvw is the number of edges in the shortest path between the nodes v and w in the graph. If v and w are in dierent connected components of the graph the distance between them is set to innite. In a directed graph it is typical to consider the directed distance ⃗ dvw between a pair of nodes v and w as the length of the directed shortest path from v to w. However, in general ⃗ dvw ̸ = ⃗ dwv, which violates the symmetry property of a metric, so that ⃗ dvw is not a distance but a pseudo-distance or a pseudo-metric.

Denition 7 A (directed) graph is (strongly) connected if there is a (directed) path between each pair of nodes. The (strongly) connected components of a (directed) graph are its maximal (strongly) connected subgraphs.

Denition 8

The following are denitions of several kinds of graphs used in this

Chapter:

An undirected tree of n nodes is a graph which is connected and has no cycles. The simplest tree is the path Pn. The path (also know as linear path or chain) is the tree of n nodes, n -2 of which have degree 2 and two nodes have degree A graph is bipartite if its nodes can be split into two disjoint (non-empty) subsets

V 1 and V 2 and V 1 ∪ V 2 = V , such that if (v, w) ∈ E, then v ∈ V i and w ∈ V j for i ̸ = j.
Bipartite graphs do not contain cycles of odd length.

A graph is planar if it can be drawn in a plane in such a way that no two edges intersect except at a node with which they are both incident.

Erd®s-Rényi (ER) G (n, p) [START_REF] Erdos | On the evolution of random graphs[END_REF] graph with n nodes: constructed by connecting nodes randomly in such a way that each edge is included in G (n, p) with probability p independent from every other edge.

Barabási and Albert (BA) one [START_REF] Barabási | Emergence of scaling in random networks[END_REF]: created on the basis of a preferential attachment process. The graph is constructed from an initial seed of m 0 vertices connected randomly like in an Erd®s-Rényi G (n, p). Then, new nodes are added to the network in such a way that each new node is connected to c ≤ m 0 of the existing ones with a probability that is proportional to the degree of these existing nodes.

Conguration model (see Chapter 7 in [237]): with a given degree sequence of length n, create n vertices and assign a degree to each vertex. Create kv half-edges or stubs at the node v according to the degree kv assigned to it. The sum of such stubs must be even. Connect pairs of stubs randomly until all of them have been used. The result is a multi-pseudograph, i.e., it may contain multiple edges and self-loops. However, for extremely large graphs the expected number of self-loops and multi-edges goes to zero.

Denition 9 Let M ∈ C n×n be a matrix. The spectrum of M is the set of not necessarily dierent numbers κ j , j = 1, . . . , n such that

M ϑ j = κ j ϑ j , (2)
where the vectors ϑ j are known as eigenvectors of M and the scalar κ j are the corresponding eigenvalues.

The following matrices and their spectra will be considered (Table 1):

name symbol denition spectrum adjacency A A ij = 1 (i, j) ∈ E 0 (i, j) / ∈ E λ 1 ≥ • • • ≥ λn Laplacian L L ij =    -1 (i, j) ∈ E k i i = j 0 otherwise 0 = µ 1 ≤ • • • ≤ µn distance D D ij = d ij i ̸ = j 0 i = j σ 1 ≥ • • • ≥ σn
Table 1: Denition of some matrices used in this paper.

5 Viral proteins 5.1 Why are proteins important for viruses?

Viruses are complex replication machines. A great part of their success relies in the functions of viral proteins. A group of viral proteins, known as structural proteins

[253], form the viral capside [266], which represents a protective coat around the viral genome, and dene the shape of the virus which can be helical, icosahedral, or complex [200]. Some viruses also contain glycoproteins at their surface, which are responsible for the attachment and penetration of the virus into the target cells [START_REF] Banerjee | Viral glycoproteins: biological role and application in diagnosis[END_REF], and consequently are vital for viral infection. In viruses containing lipid membranes at their capside, the matrix proteins are needed by the virus after it has entered a host cell to expel its genetic material. The second legion of viral proteins are found inside the virus capside and are known as nonstructural (NS) proteins [253]. These proteins develop their functions inside the infected cell where they help in virus replication and virus assembly. Other NS proteins play their role as immunomodulators, that is by minimizing the immune system's targeting of the virus. The viral proteins army is completed with accessory and regulatory proteins [253]. They play diverse roles in several viral functions, ranging from transcription rate regulation of viral genes encoding structural proteins to modication of host cell functions.

For these reasons, all viruses contain proteins, ranging from only a few kinds like in the case of human papilloma virus which has 6-8 proteins, to up to approximately 1000 proteins contained by other viruses like the Acanthamoeba polyphaga mimivirus. All viruses are obligate intracellular parasites, which means that their replication requires to sequester host cell functions [241,242]. During their sequestration of the host cell machinery viruses also neutralize innate host defenses that are aimed at inhibiting the capacity of infected cells for protein synthesis. These functions result in a complex system of interactions between virus and host proteins making that almost every step of the translation process is targeted by virally encoded functions.

General denitions

We start by introducing some general concepts about proteins and their structures

[250].

Denition 10 A protein P is a polymer in which the monomeric units are amino acids, A . An amino acid is an organic compound having the groups -NH2 (amino) and COOH (carboxylic acid) in the same structure. For proteins found in nature both groups are bonded to the same carbon atom, which is denominated as alpha-carbon,

Cα.

A third chemical group is attached to this Cα, which we represent generically as R. The sequence of amino acids in a protein denes its primary structure. The list of amino acids and their 3-and 1-letter codes are given in Table 2.

Denition 11 Let P be a protein with A amino acids. The backbone of P is the linear chain formed by the sequence (N -Cα -C) A , where N and C are the nitrogen and carbon atoms of the amino and carbonyl groups, respectively, and Cα is the carbon atom supporting the dierent substituents attached to each dierent amino acid. The three dimensional structure of the backbone of a protein denes its secondary structure.

It is well-known that the bond lengths between a pair of atoms in a molecule is almost constant around an average value with a small deviation due to the thermal eects.

The same happens for the bond angles formed by triads of atoms. Therefore, the largest variability in a protein backbone is given by the dihedral angles formed between tetrads of atoms.

P (A 1 , A 2 , A 3 ) and Q (A 2 , A 3 , A 4 ).
Obviously, 0 • ≤ φ ≤ 180 • , where the lower dihedral angle occurs when the four atoms are folded in the same plane, a rotation of the planes arrives at φ = 90 • where the two planes are perpendicular, until it arrives at the value of φ = 180 • where the four atoms are in an unfolded state as illustrated in Fig. 10. Let us now dene an important concept for the denition of a folding degree index.

Denition 13 Two dihedral angles φ 1 and φ 2 are adjacent if they share exactly one plane between them.

An example is provided in Fig. 11 where the angles φ 1 and φ 2 are adjacent because they share the plane P 2 between them.

We now consider a linear chain of adjacent dihedral angles [START_REF] Estrada | Characterization of 3D molecular structure[END_REF][START_REF] Estrada | Characterization of the folding degree of proteins[END_REF].

Fig. 11: Illustration of the adjacency between two dihedral angles. The rst dihedral angle φ 1 is formed between the planes P 1 and P 2 and the second φ 2 is formed between the planes P 2 and P 3 .

Denition 14 Let P be a protein with A amino acids. The graph of adjacent dihedrals (GAD) of P is the path graph Pn = (V, E) with n = 3A -3 where the ith vertex in V corresponds to the ith dihedral angle in the backbone of P and two nodes i and j are connected, i.e., (i, j) ∈ E, if and only if the corresponding dihedral angles are adjacent.

Remark 2 The structural information needed to built the dihedral angles of the protein backbone is reported on the Protein Data Bank (PDB) [START_REF] Berman | The protein data bank[END_REF]. In the PDB the structural information about every atom is stored in the form illustrated in Fig. 12, where the third column indicates the name of the atom in the corresponding amino acid with name provided in the fourth column and number in the sixth one, of the protein chain indicated in the fth column. For instance, N GLY A 1 indicates that this entry corresponds to the nitrogen atom of the glycine 1 of the chain A. The columns 7-9 give the x, y, z coordinates of this atom in the crystal as determined by x-rays crystallography.

Main protease of coronaviruses

The 29 proteins produced by SARS CoV-2 are described in Table 3 (see [100]). They constitute potential pharmacological targets for drug repurposing. Here we will focus on the main protease and compare it with that of SARS-CoV-1. Viral proteases are enzymes that catalyze the maturation processing of polyproteins in the host cell interior such that their catalytic activity is required for the production of new, infectious virions [108,232]. Therefore, they are attractive pharmacological targets for inhibiting viruses of diverse nature as it was demonstrated for instance for HIV protease inhibitors [232].

The main protease of coronaviruses plays a fundamental role in virus replication as it is responsible for the proteolytic cleaves of the overlapping polyproteins to functional proteins. The main protease of SARS-CoV-1 was already proposed as a target for drugs inhibiting the coronavirus producing the epidemics of 2002-2003 [12, 255]. Therefore it is interesting to investigate how dierent the main proteases of SARS-CoV-1 and SARS-CoV-2 are towards the potential investigation of drugs that act across dierent variations of coronaviruses. We consider here the three-dimensional structures of these proteases from the PDB [START_REF] Berman | The protein data bank[END_REF]. For the case of M pro of SARS CoV-1 we study here the structures with PDB code: 4).

Before starting with our analysis let us explore some of the classical ways of comparing the structures of the two main proteases. First, we consider the amino acid sequences of both proteases (here done for 2BX4 and 6Y2E for the sake of simplicity).

As can be seen in Fig. 13 there are change in only 12 amino acids out of 306, which represents 96% of similarity between the two sequences.

Additionally, we perform a superposition of the 3-dimensional structures of both proteases as illustrated in Fig. 14. The root mean square distance (RMSD) between the coordinates of the the two structures is only 0.72Å, which indicates that both structures display high 3D similarities.

In the following we will investigate whether networks analysis is able to nd any signicant dierences in the structures of the two main proteases which may point out towards some evolutionary hints separating both coronaviruses.

6.1 Folding degree 6.1.1 What is folding?

Let us start with a simple example. We consider a linear chain formed by seven straight segments of the same length. The angle between two segments is xed at 90 previously used for modeling a linear chain that shows some protein-like properties [147].

In Fig. 15 (b) and (c) we have two embeddings of the same linear chain into a cube by folding the planes as indicated before. We can say that these two linear chains are folded as they are not in the totally extended form given in Fig. 15 (a). The structure in Fig. 15 (b) has three dihedral which fully-folded, e.g., which has 0 • , and two others have 90 • . Contrastingly, in Fig. 15 (c) there is only one dihedral which is fully-folded while the rest of pairs of adjacent planes are perpendicular. Then, it is obvious that the folding degree of the two linear chains is necessarily dierent, with structure (d) being more folded than structure (e). Last but not least, the vertices of the two linear chains (intersections of adjacent segments) occupy the 8 vertices of the cube. This means that the sum of the distances between the center of the cube to each vertex of the linear chain are identical for both chains in Fig. 15. In other words, the radius of gyration of both structures, which account for the their respective packing, are exactly the same. Therefore, there are two important conclusions that can be extracted from this example:

1. The folding degree of linear chains where the segment length and inter-segment angles are xed depends only on the dihedral angles between the planes formed by adjacent pairs of segments;

2. The folding degree of a linear chain does not depend on the packing of the corresponding structure.

Folding degree of proteins

The folding degree index was proposed and studied by Estrada [6366] and it has been recently extended by Sladek et al. [217]. Let Pn be the DAG of a protein P. Let i and j be the ith and jth dihedral angles in Pn. We then consider a particle which is created at the dihedral i after being annihilated at the dihedral j. Let c † i (c j ) be the creation (annihilation) operators of this particle at the corresponding dihedrals and let t ij be the energy needed to transfer such particle between the two dihedrals [START_REF] Estrada | Tight-binding 'dihedral orbitals' approach to electronic communicability in macromolecular chains[END_REF]224].

Then, the energy of the DAG can be determined by using the following tight-binding Hamiltonian [101]: H = -

(i,j)∈E t ij c † i c j + t ji c † j c i + i∈V V i c † i c i , (3)
where we have considered that the transfer of the particle may occur in both directions and where V i is the potential energy needed to move the particle to the corresponding dihedral.

Our rst assumption is to have t ij = t ji = -1, where the negative sign indicates that it cost energy to move a particle from one dihedral to another. Let us now think about V i . Obviously, it has to be a property of the corresponding dihedral angle. When the dihedral angle is zero it costs very little to move the particle from one dihedral to the other (see Fig. 10). We could expect a smooth increase of this cost from φ j = 0 • to φ j = 90 • due to the increase in the separation between the two planes. However, when the angle is 90 • < φ j ≤ 180 • we consider that the two planes forming the dihedral are in an extended or unfolded conformation, such as the potential function is set to its maximum. To capture this behavior we set [6366]

V i = - 1 2 [1 + sgn (cos (φ i ))] cos (φ i ) , (4)
as the potential function for the dihedral angle, which is illustrated in Fig. 16.

Consequently, the Hamiltonian has the following appearance: which immediately implies that H = -W , where W is a vertex-weighted adjacency matrix of the GAD.

H =       -V 1 -1 • • • 0 -1 -V 2 . . . . . . . . . . . . . . . -1 0 -1 -V k       , (5)
The Boltzmann probability of nding the system in a state with energy E j when the inverse temperature of the system is τ = (k B T ) -1 with k B being a constant and T is the temperatureτ is typically represented by β in statistical physics, but this letter is already reserved here for a dierent variable, is

p j (τ ) = e -τ Ej Z , (6)
where Z = tr e -τ H . Therefore, the Boltzmann probability of the system is given by

p j (τ ) = e τ λj I 3 (G, τ ) , (7)
where

I 3 (G, τ ) := Tre τ W . ( 8
)
From the own denition of the index I 3 we can see that it plays the role of the partition function of the GAD. We now can dene the entropy of the graph as [START_REF] Estrada | Statistical-mechanical approach to subgraph centrality in complex networks[END_REF] S (G, From the general expression of the entropy one can obtain the graph enthalpy H (G, τ ) =j p j λ j and the graph free energy, which is sometimes named the natural connectivity of the network [START_REF] Estrada | Statistical-mechanical approach to subgraph centrality in complex networks[END_REF]:

τ ) = -k B j p j ln p j = - 1 T j (p j λ j ) + k B ln I 3 (G, τ ) .
F (G, τ ) = -τ -1 ln I 3 (G, τ ) . (10)
Therefore, when τ is taken to be unit, the free energy of the GAD is given by minus the logarithm of the I 3 index.

As an illustration we consider the two linear chains shown in Fig. 15 (c) and (d) and calculate the corresponding I 3 indices. For the structure (c) we have I 3 ≈ 22.753 and for (d) it is I 3 ≈ 13.378. Thus, as expected the structure (c) of Fig. 15 is more folded than (e). As a second example we consider here the same protein, the ribonuclease-A, which was studied at nine dierent temperatures from 98 to 320 K [230]. It is known that increasing the temperature produce that the 3-dimensional structure of the protein unfolds, giving rise to its denaturation. However, such changes in the range of temperatures studied is not too large as can be seen in Fig. 17 On the structural meaning of I 3 (G, τ )

Here we provide a structural interpretation of the folding degree index which have not been previously discussed in the literature. The index I 3 (G, τ ) can be expressed as:

I 3 (G, τ ) = ∞ k=0 Tr τ k W k k! . (11)
Let us write down the expression for some of the spectral moments of the matrix W :

Tr (W ) = i cos φ i (12) Tr W 2 = 2m + i cos 2 φ i (13) Tr W 3 = 3 i k i cos φ i + i cos 3 φ i (14)
Tr

W 4 = 2m + 4P 3 + 4 (i,j)∈E cos φ i cos φ j + 4 i k i cos 2 φ i + i cos 4 φ i (15) Tr W 5 = 5 i (5k i + 4δ i + 4ι i ) cos φ i + 5 i k i cos 3 φ i + 5 (i,j)∈E cos 2 φ i cos φ j + cos φ i cos 2 φ j + i cos 5 φ i , ( 16
)
where m is the number of pairs of adjacent dihedrals in the GAD, k i is the number of adjacent dihedrals that i has, i.e., either 1 for the terminal angles or 2 for the rest, P 3 is the number of triples of adjacent dihedrals, δ i =1 terminal angles and the ones immediately adjacent to them or δ i =2 for the rest, and ι i = 0 for terminal dihedrals or ι i = 1 for the rest.

As can be seen these expressions contain information not only about the individual dihedral angles, like in i cos φ i but also about pairs of adjacent dihedrals like in (i,j)∈E cos φ i cos φ j , and other structures which appear in higher order moments like pairs of dihedrals separated by a third one which appear in Tr W 6 , triples of adjacent dihedrals which appear in Tr W 7 , and so on.

All of these contributions are due to the fact that the index I 3 (G) counts the number of weighted closed walks in the GAD. The connection between walks and the powers of the adjacency matrix comes from the following well-known result.

Theorem 1 Let G be a simple graph with adjacency matrix A. Let v, w ∈ E, then the number of walks of length k between the nodes v and w is given by A k vw .

Then, every closed walk of a given length encloses a weighted linear connected subgraph of the GAD as shown partially by the expressions of the rst moments of the adjacency matrix of the GAD shown before. The inclusion of this nonlocal information about dihedral angles is revealed in the following example. In Fig. 18 we illustrate the structures of six dierent conformations of a linear chain of six atoms. Therefore, there are n = 3 dihedral angles in this chain such that we can form GADs of three nodes and the corresponding cos φ i in the main diagonal. As can be seen in Table 5 the sum of the dihedral angles is not enough to dierentiate these conformers, which are very well characterized by the index I 3 (G) . Indeed, this index increases its value with the increase of the folding degree of the linear chain, therefore it is a measure of the degree of folding of such chains. Table 5: Values of the sum of the cosines of the dihedral angles and folding degree indices of the conformations of a linear chain illustrated in Fig. 18.

Mathematical properties of the folding degree index

Let P k be a path graph with k nodes. Let A (P k ) be the adjacency matrix of P k . Then, if a protein has N amino acids we have that W = A (Pn) + V , where n = 3N -3, V is a diagonal matrix having the values of V i in its main diagonal for the corresponding dihedral angles.

Before proceeding we introduce the following. Denition 16 Let α (M ) be the condition number of M . Then, the Specht ratio at α (M ) is [34] S (α (M )

) := (α (M ) -1) α (M ) 1/(α(M )-1) e log α (M ) . ( 17
)
Let us rst prove the following auxiliary result.

Lemma 1 Let α e A(Pn) and α e V be the condition number of the matrices exp (A (Pn)) and exp (V ), respectively. Then S α e A(Pn) > S α e V , [START_REF] Banerjee | Viral glycoproteins: biological role and application in diagnosis[END_REF] for any linear protein with N ≥ 2 residues.

Proof In any protein with N ≥ 2 residues, n = 3N -3 ≥ 3. The condition number of

exp (V ) is α e V = max i e λi(V ) min i e λi(V ) = max i e cos φi min i e cos φi ≤ e 1 e -1 = e 2 . ( 19
)
In the case of exp (A (Pn)) we have α e A(Pn) = max i e λi(A)

min i e λi(A) = exp 2 cos π n+1 exp -2 cos π n+1 = exp 4 cos π n+1 ≥ exp √ 8 > e 2 . ( 20
)
Then because the Specht function is increasingly monotone we have the nal result [START_REF] Bourin | Reverse inequality to GoldenThompson type inequalities: Comparison of ea+ b and eaeb[END_REF].

Theorem 2 Let α (V ) be the condition number of V . The folding degree index is bounded as

e λn(V ) S (α (V )) EE (Pn) e λn(A(Pn)) S (α (V )) Tre V        ≤ I 3 (G) ≤ e λ1(V ) EE (Pn) e λ1(A(Pn)) Tre V . ( 21
)
Proof Let us rst prove the upper bounds. Using the well-know Golden-Thompson (see [START_REF] Bushell | Trace inequalities for positive denite matrix power products[END_REF][START_REF] Fang | Inequalities for the trace of matrix product[END_REF]) inequality we have

I 3 (G) = Tre A(Pn)+V ≤ Tr e A(Pn) e V . (22)
Then, because both matrices exp (A (Pn)) and exp (V ) are positive denite we have Tr e A(Pn) e V ≤ e λ1(V ) Tre A(Pn) e λ1(A(Pn)) Tre V .

(

) 23
For the lower bounds, we have that Pn) .

I 3 (G) = Tre A(Pn)+V ≥ S -1 (α (V )) Tr e V e A(
(

) 24
We then use the fact that Tr e V e A(Pn) ≥ e λn(V ) EE (Pn)

e λn(A(Pn)) Tre V , (25)
to prove the nal result.

⊓ ⊔

Corollary 1 Let n → ∞, then the index I 3 (G) is bounded as

(n -1) -2 cosh (2) e 3 Tre V e 4      ≤ I 3 (G) ≤ ne -e -2e cosh (2) e 2 Tre V . ( 26
)
Corollary 2 Let n → ∞, and let φ = 1 n i φ i . Then

I 3 (G) ≥ ne φ-2 e 2 . ( 27
)
Lemma 2 The folding degree index is bounded as

I 3 (G) ≤ j e cos φj +2 cos    πj n + 1    . ( 28
)
Proof The result follows from the Golden-Thompson inequality and the fact that

Tr e V e A(Pn) ≤ j e λj (V ) e λj (A(P n )) . (29)
⊓ ⊔

Folding degree of SARS-CoV main proteases

Although the length of the proteases is 306 amino acids there are structures which are only resolved for amino acids 3 to 300, which gives a length of 298 (see Table 4). Thus, for the sake of homogeneity of the analysis we consider hereafter the same part of the amino acids sequence for all the structures analyzed, i.e., from residue 3 to residue 300.

This does not alter the analysis as the two extremes of the protease are disordered and do not participate in important interactions. In Table 6 we give the values of the folding degree index I 3 (G) /N , the folding degree entropy S (G) and free energy F (G). As can be seen the mean of the folding degree of both classes of proteases are very close to each other. Indeed, a two-sample F-and t-test for equal variances conrm that both samples comes from normal distributions with the same variance. The probability that they are dierent is only 0.69. For the cases of the folding entropy and free energy the similarities are even bigger, with probabilities that they come from dierent normal distributions of only 0.203 and 0.212, respectively. Then, clearly, the degree of folding of the main protease of SARS-CoV-2 is not dierentiable from that of SARS-CoV-1, indicating that any evolutive processes from SARS-CoV-1 to SARS-CoV-2 has not changed this parameter in a signicant way. Table 6: Values of the folding degree index, the folding degree entropy and free energy of the main protease of SARS-CoV-1 and SARS-CoV-2 studied here.

SARS CoV-1 SARS CoV-2 PDB I 3 (G) /N S (G) F (G) PDB I 3 (G) /N S (G) F (G) 2H2Z
Obviously, it is possible that the global parameter used before hides some important changes in the main protease occurring at a local level. For that reason we investigate here the amino acid contributions to the global folding degree of the main proteases of SARS-CoV-1 and SARS-CoV-2. For that we dene the contribution of the ith dihedral angle of a given protein to the global I 3 (G, τ ) index as

I 3 (i, τ ) := e τ W ii . ( 30
)
Notice that I 3 (G, τ ) = N -3 i=1 I 3 (i, τ ). Then, we take the average of the contribu- tions made by the dihedral angles forming a given amino acid as its contribution to the global I 3 (G, τ ) index [65]. In Fig. 19 we illustrate the contributions of the amino acids in the chains of the main proteases studied here for SARS-CoV-1 (a) and SARS-CoV-2 (b). In general, there is a good agreement in the amino acid contributions for each of the proteases in the same class, i.e., SARS-CoV-1 and SARS-CoV-2. As our goal is to see if there are signicant dierences between the two classes of proteases we use the sequences of amino acid contribution to the folding degree as a vector to classify the proteins into dierent structural classes using a binary hierarchical cluster. The results are illustrated in Fig. 20. It can be seen that there is a great similarity between the structures 6M2Q and 6Y2E which correspond to SARS-CoV-2.

However, the three SARS-CoV-2 proteases are not signicantly dierent from the pair 1UJ1 and 2H2Z, which correspond to SARS-CoV-1. Therefore, in general the amino acid contributions to the folding degree does not reect any signicant change from SARS-CoV-1 and SARS-CoV-2, indicating that neither the global nor the local folding degree of the protease has been a signicant target of the evolution of this protein. of SARS-CoV-2, using a binary hierarchical cluster.

Protein residue networks

Another way of representing a protein is through the interaction between its amino acids (residues). These networks are know as protein residue networks (PRN) (see [START_REF] Estrada | The structure of complex networks: theory and applications[END_REF] Chapter 14, [122, 214]). The nodes of the PRN represent the α-carbon of the amino acids, which are the ones labeled as CA in the PDB le. Then, we consider a cuto radius r C , which represents an upper limit for the separation between two residues in contact. The distance r ij between two residues i and j is measured by taking the distance between Cα atoms of both residues. Then, we have the following.

A ij = H (r C -r ij ) i ̸ = j, 0 i = j, (31)
where H (x) is the Heaviside function. This overall three-dimensional arrangement of a protein chain in space denes its tertiary structure.

Here we use r C as the typical interaction distance between two amino acids, which is equal to 7.0 Å. We have tested distances below and over this threshold obtaining in general networks which are either too sparse or too dense, respectively. Denition 18 Let G be a simple graph. A subgraph ngerprint of G is a vector

F = (F 1 , F 2 , • • • , F k )
where F i counts the number of subgraphs of a given type in G.

Here we consider a subgraph ngerprint composed by the 21 subgraphs illustrated in Fig. 22, whose numbers can be obtained analytically as follow (formulae for subgraphs

F i (i = 1, • • • 20
) are adapted from [START_REF] Alon | Finding and counting given length cycles[END_REF] and for F 21 has been obtained here).

Theorem 3 Let k i and t i be the degree and the number of triangles at the node i.

Then, the number of subgraphs F i illustrated in Fig. 22 are obtained as follow:

F 1 = n, (32)
F 2 = 1 2 i k i , (33)
F 3 = 1 2 i k i (k i -1) , (34)
F 4 = 1 6 tr A 3 , (35)
F 5 = (i,j)∈E (k i -1) (k j -1) -3F 4 , (36)
F 6 = 1 6 ki≥3 k i (k i -1) (k i -2) , (37)
F 7 = 1 8 tr A 4 -4F 3 -2F 2 , (38)
F 8 = ki>2 t i (k i -2) , (39)
F 9 = 1 4 i,j A 2 ij A ij A 2 ij • A ij -1 , (40)
F 10 = 1 10 tr A 5 -30F 4 -10F 8 , (41)
F 11 = 1 2 ki≥4 t i (k i -2) (k i -3) , (42)
F 12 = 1 2 ki>2 (k i -2) × i,j A 2 ij 2 -2F 9 , (43)
F 13 = (i,j)∈E A 2 ij (k i -2) (k j -2) -2F 9 , (44)
F 14 = i t i   i̸ =j (A 2 ) ij   -6F 4 -2F 8 -4F 9 , (45)
F 15 = 1 2 i t i (t i -1) -2F 9 , (46)
F 16 = (i,j)∈E A 3 ij A 2 ij -9F 4 -2F 8 -4F 9 , (47)
F 17 = 1 12 tr A 6 -2m -12F 3 -24F 4 -6F 5 -12F 6 -48F 7 -36F 9 -12F 12 -24F 15 , (48)
F 18 = 1 2 ki>2 (k i -2) B i -2F 16 , (49)
where

B i = A 5 ii -20t i -8t i (k i -2)-2 (i,j)∈E A 2 ij (k j -2)-2 (i,j)∈E t j -A 2 ij , (50)
F 19 = (i,j)∈E A 2 ij 3 , (51)
F 20 = i t i • i̸ =j A 2 ij 2 -6F 9 -2F 16 -6F 19 , (52)
F 21 = 1 14 tr A 7 -126F 4 -84F 8 -112F 9 -70F 10 -28F 11 -14F 13 -14F 14 -56F 16 -14F 18 -84F 19 -28F 20 ) . (53)

Network ngerprints of SARS-CoV main proteases

Using vectors of network ngerprints we study again the classication of the seven structures of the main protease of SARS-CoV-1 and SARS-CoV-2 with a binary hierarchical cluster. The results illustrated in Fig. 23 reveal a high similarity between the structures 2BX4 and 1UJ1, both corresponding to SARS-CoV-1 as well as between 2H2Z and 2DUC (SARS-CoV-1) and 6M2Q and 6M03 (SARS-CoV-2). However, the structure 6Y2E which belongs to SARS-CoV-2 is more similar to SARS-CoV-1 proteases than to SARS-CoV-2 ones. All in all there is not a clear indication about a global signicant dierentiation between the two classes of proteases which points out to an evolutionary route from SARS-CoV-1 to SARS-CoV-2.

However, when we analyze the relative dierence between the individual subgraphs in both groups of proteases we observe a very interesting pattern. In Fig. 24 we illustrate these results where the label of the subgraphs is the one given in Fig. 22. It can be This could be an indication of the existence of more potential cavities in SARS-CoV-2 than in its predecessor. We will go back to this important problem later in this chapter.

Global network descriptors

Here we consider a series of classical global network structure descriptors for the analysis of the main proteases of SARS-CoV-1 and SARS-CoV-2. Complete description,

properties and references about these descriptors are found in [START_REF] Estrada | The structure of complex networks: theory and applications[END_REF]. The rst category of measures correspond to those related to the most local structure around the nodes, such as those based on the degree of the nodes. The degree accounts for the immediate eect of a node to its closest neighborhood. Among these measures we use here the edge density, which is dened as Denition 19 Let G be a simple graph with n vertices and F 2 edges. Then, the edge Lemma 3 Let G be a simple connected graph with n vertices. Then, the edge density of G is bounded as

density of G is δ (G) = F 2 ( n 2 ) = 2F 2 n (n -1)
2 n = δ (Tn) ≤ δ (G) ≤ δ (Kn) = 1, ( 55
)
where Kn is the complete graph with n nodes and Tn is any tree with n nodes.

Remark 3 If non-connected graphs are allowed, the lower bound of δ (G) is zero, which is reached for the trivial graph Kn.

Denition 20 Let G be a simple connected graph with n vertices and let d (i, j) be the shortest path distance between the vertices i and j. Then, the average shortest path distance of G is

L (G) = 1 n (n -1) i<j d (i, j) . ( 56
)
Lemma 4 Let G be a simple connected graph with n vertices. Then, the average shortest path distance of G is bounded as

1 ≤ L (G) ≤ n -1 3 , ( 57
C (G) = 3F 4 F 3 . ( 58
) Remark 4 Let M 1 (G) = i k 2
i be the rst Zagreb index (see [263] and refs. therein) of G. Then,

C (G) = F 4 M 1 (G) -2F 2 . ( 59
)
Denition 22 Let G be a simple graph and let t i be the number of triangles attached to the node i.

C (G) = 1 n n i=1 2t i k i (k i -1) . ( 60
)
Remark 5 Both indices C (G) and C (G) are bounded between zero and one, where the lower bound is obtained for any triangle-free graph and the upper one is attained for any maximally clustered graph. For important dierences between them, the reader is referred to [START_REF] Estrada | When local and global clustering of networks diverge[END_REF]. Denition 23 Let G be a simple connected graph with n vertices and let k i be the degree of the node i. Then, the degree assortativity of G is the Pearson correlation coecient of the degree-degree correlation of G

[177] r (G) = 1 m (i,j)∈E k i k j - 1 2m (i,j)∈E (k i + k j ) 2 1 2m (i,j)∈E k 2 i + k 2 j - 1 2m (i,j)∈E (k i + k j )
2 .

(

) 61
Remark 6 r > 0 indicates degree assortativity and r < 0 indicates degree disassortativity.

Theorem 4 [START_REF] Estrada | Combinatorial study of degree assortativity in networks[END_REF] Let G be a simple connected graph. Then, the degree assortativity coecient of G is

r (G) = F 3 F 5 F 3 + C (G) - F 3 F 2 3F 6 + F 3 1 - F 3 F 2 . ( 62
)
Denition 24 Let G be a simple connected graph. The degree heterogeneity index of G is dened as [START_REF] Estrada | Quantifying network heterogeneity[END_REF][START_REF] Estrada | Degree heterogeneity of graphs and networks. i. interpretation and the "heterogeneity paradox[END_REF][START_REF] Estrada | Degree heterogeneity of graphs and networks. ii. comparison with other indices[END_REF] ρ (G) =

(i,j)∈E k -1/2 i -k -1/2 j 2 . ( 63
)
Theorem 5 [START_REF] Estrada | Quantifying network heterogeneity[END_REF][START_REF] Estrada | Degree heterogeneity of graphs and networks. i. interpretation and the "heterogeneity paradox[END_REF] Let G be a simple connected graph with n nodes. The degree heterogeneity index of G is bounded as

0 ≤ ρ (G) ≤ n - √ n -1, (64)
where the lower bound is attained for any regular graph and the upper bound is reached if and only if G ∼ = Sn, where Sn is the star graph with n nodes.

Denition 25 Let G be a simple connected graph with n nodes. Let ρ ikj be the number of shortest paths between the nodes i and j that cross the node k, and let ρ ij be the total number of shortest paths that go from i to j. Then, the betweenness centrality of the node k is [START_REF] Freeman | Centrality in social networks conceptual clarication[END_REF] BC

k = i̸ =k̸ =j ρ ikj ρ ij , ( 65
)
and BC is the average betweenness of the nodes of G.

Denition 26 A walk of length k in G is a set of nodes i 1 , i 2 , . . . , i k , i k+1 such that for all 1 ≤ l ≤ k, (i l , i l+1 ) ∈ E. A closed walk is a walk for which i 1 = i k+1 .
Theorem 6 Let G be a simple graph with adjacency matrix A. Then, the number of walks of length k between the nodes i and j in a network is given by A k ij .

Denition 27 Let G be a simple graph. The subgraph centrality of the node i is dened as [START_REF] Estrada | Subgraph centrality in complex networks[END_REF] SC

i = ∞ k=0 A k ii k! = (exp (A)) ii , (66)
and SC is the average subgraph centrality of the nodes of G.

Remark 7

The average subgraph centrality of the nodes of G is SC = 1

n Tr e A where Tr e A =: EE (G) is known as the Estrada index of G [54].

Theorem 7 [START_REF] De La Peña | Estimating the Estrada index[END_REF] Let G be a simple connected graph with n nodes. Then

(n -1) -2 cosh (2) ≤ EE (G) ≤ e n-1 + (n -1) e -1 , (67)
where the upper bound is attained if and only if G ∼ = Kn and the lower one is obtained

if G ∼ = Pn when n → ∞.
Theorem 8 [START_REF] Benzi | Quadrature rule-based bounds for functions of adjacency matrices[END_REF] Let G be a simple graph and let a, b ∈ R be such that the spectrum of A is contained in [a, b]. Then, the Estrada index of G is bounded as

n i=1 b 2 exp (k i /b) + k i exp (-b) b 2 + k i ≤ EE (G) ≤ n i=1 a 2 exp (k i /a) + k i exp (-a) a 2 + k i , (68)
where k i is the degree of the node i. Denition 28 Let G be a simple graph. The communicability between the nodes i and j is Γ ij = (exp (A)) ij and Γ (G) is the average communicability between all pairs of dierent nodes in G [START_REF] Estrada | Communicability in complex networks[END_REF][START_REF] Estrada | The physics of communicability in complex networks[END_REF].

Lemma 5 Let Γ (G) be the average communicability between all pairs of dierent nodes in a simple graph G with n nodes. Then,

Γ (G) ≤ EE (G) n = SC, (69)
where EE (G) is the Estrada index of G.

Proof Let us write Γ (G) as

Γ (G) = 1 n (n -1)
Tr e A J -Tr e A , [START_REF] Estrada | The communicability distance in graphs[END_REF] where J is the corresponding all-ones matrix. Then, because both J and exp (A) are Hermitian and positive semidenite following [START_REF] Bushell | Trace inequalities for positive denite matrix power products[END_REF][START_REF] Fang | Inequalities for the trace of matrix product[END_REF] we have Tr e A J ≤ λ 1 (J) Tr e A = nEE (G) .

Γ (G) ≤ 1 n (n -1) [nEE (G) -EE (G)] = EE (G) n = SC. (72)
⊓ ⊔

Theorem 9 [START_REF] Estrada | The communicability distance in graphs[END_REF][START_REF] Estrada | Communicability angle and the spatial eciency of networks[END_REF] Let G be a simple graph and let:

ξ ij := SC i + SC j -2G ij , (73)
θ ij = cos -1 G ij G ii G jj . (74)
Then, ξ ij is a Euclidean distance between the vertices i and j, and θ ij is the Euclidean angle spanned between the position vectors of the nodes i and j embedded in a hypersphere. Remark 9 We consider here the averages ξ and θ of the distances and angles between all pairs of dierent nodes in a graph G.

Denition 29 Let G be a simple connected graph with adjacency matrix A. Let λ 1 > λ 2 ≥ • • • ≥ λn be the eigenvalues of A. Then, the ith entry of the vector ψ 1 , such that Aψ 1 = λ 1 ψ 1 , is the eigenvector centrality [START_REF] Bonacich | Power and centrality: A family of measures[END_REF] of the node i and ĒC is the average eigenvector centrality of the nodes of G.

Theorem 10 [START_REF] Cvetkovic | Eigenspaces of graphs[END_REF] Let G be a non-bipartite graph. Let N k (i) be the number of walks of length k that start at the node i. Then,

EC i = lim k→∞ N k (i) n j=1 N k (j) . ( 75
)
Denition 30 Let G be a simple connected graph. Let G i ⊂ G be a cluster of nodes and edges of G, such as δ (G i ) ≫ δ ( Ḡi ). Then, G i is said to be a community in G.

Then, for a given partition of G into n C communities, the modularity index Q is dened

as [178] Q = n C k=1   |E k | m - 1 4m 2   n j=1 k j   2   , (76)
where |E k | is the number of edges between nodes in the kth community of the network, m is the total number of edges in the network and k j is the degree of the node j.

Global descriptors of SARS-CoV main proteases

Here we have calculated the 13 descriptors previously dened for the 7 PRNs of the main proteases of SARS-CoV-1 and SARS-CoV-2. In Table 7 we report the average and standard deviations of the four SARS-CoV-1 proteases as well as of the three SARS-CoV-2 ones. Then, we report the relative dierences between these means for each descriptor and applied the two samples F -test and Kolmogorov-Smirnov (KS)

test. The values reported in Table 7 for these tests correspond to the probabilities that the two samples come from populations with dierent variances. As can be seen, none of the descriptors dierentiate the two samples according to the F -test (p > 0.95) and only the descriptors based on exp (A) dierentiate the two samples according to the KS test. In order to interpret these results and prepare the land for the next steps let us state the following result.

Lemma 6 Let G be a simple graph. Then, the Estrada index of G is bounded as

EE (G) ≥ F
F 21 . ( 77
)
Proof Based on the relations shown before for tr A k for k ≤ 7 we have that the right-hand-side part of eq. ( 77) is

7 k=0 tr(A k ) k!
from which the inequality follows. ⊓ ⊔ Accordingly, it is obvious that the descriptors based on exp (A) count the number of subgraphs of any kind existing in the network, giving more weight to the smaller than the longer ones. For instance, if we make a truncation of the power series up to k ≤ 7 each triangle receives a weight equal to 157 126

, a square received one of 2 5

, a pentagon of 7 72 , an hexagon of 1 60

and an heptagon of 1 360

. As we have seen in the analysis of subgraph ngerprints, there are indications about an important role played by longer cycles in the main protease of SARS-CoV-2 relative to that of SARS-CoV-1.

The current results based on SC and Γ indicates that as average the nodes of the PRN of the main protease of SARS-CoV-2 participate in larger number of subgraphs than those of SARS-CoV-1. In the next subsection we investigate whether a relaxation in the penalization of longer subgraphs increases this dierentiation between both proteases.

As can be seen the structures III and IV are almost indistinguishable by Tr (exp (A)).

That is, the matrix exponential hardly dierentiate between cycles of length 5 and 7

(graph (III)) from cycles of length 6 (graph (IV)). This is due to the factorial penalization used in the power series of this matrix function, which penalizes very heavily closed walks of relatively long length [START_REF] Estrada | Generalized walks-based centrality measures for complex biological networks[END_REF][START_REF] Estrada | Accounting for the role of long walks on networks via a new matrix function[END_REF].

The Tr (exp (A)) of a graph can be written as

Tr (exp (A)) = tr ∞ k=0 A k Γ (k + 1) , (79)
due to the fact that Γ (k + 1) = k!, where Γ (•) is the Euler gamma function:

Γ (z) = x z-1 e -x dx, R (z) > 0.

(

) 80
Therefore we can generalize the subgraph centrality to E α,β (G) ii and the com- municability function to E α,β (G) ij where [1, 15, 67, 77, 85]

E α,β (G) = ∞ k=0 A k Γ (αk + β) . ( 81
)
Obviously, E 1,1 (G) = exp (A). The summation of the powers of the adjacency matrix divided by the corresponding gamma functions is the denition of the Mittag-Leer matrix functions [START_REF] Garrappa | Computing the matrix mittag-Leer function with applications to fractional calculus[END_REF] of the adjacency matrix. Below we give the relative errors between the four dierent structures illustrated in Fig. 25 for the Tr (E α,1 (G)) for α = 0.75 and α = 0.5, respectively

RD (%) = I II III IV     0.
Analysis of coronaviruses main proteaseIt is clear that as α is smaller than one the longer walks are less penalized and make more contributions to the corresponding indices. We then apply the E α,1 (G) matrix function and obtain Tr (E α,1 ) as well as Ῡ the sum of the Mittag-Leer communicability for all pair of dierent nodesin the PRN of the main proteases of SARS-CoV-1 and SARS-CoV-2 studied. The results are given in Table 8 for α = 0.75 and α = 0.5. The statistical analysis of the result is very revealing. As can be seen in Table 9 when α = 0.75 the relative dierence between the Mittag-Leer subgraph centrality and communicability of SARS-CoV-2 is more than 47% higher than those of SARS-CoV-1.

The Kolmogorov-Smirnov test indicates that the two samples come from populations with signicantly dierent variances. However, the F -test does not dierentiate between the two samples. When α = 0. Conclusions In closing, we have shown here how networks virology techniques are able to reveal some structural dierences between the main proteases of SARS-CoV-1 and SARS-CoV-2 which are not easy to nd by using traditional techniques of structural biology, such as those based on 3-dimensional superposition of structures.

Perturbation transmission as a (reaction) dynamics on networks

We give here a dynamics interpretation of the results obtained in the previous subsection. Let us consider that the nodes of a PRN G can be in either of two states: either unperturbed or perturbed by an external factor, e.g., solvent interaction, another protein, a ligand or an inhibitor. A perturbed node can transmit the perturbation to any other node in the PRN to which it interacts with. Then, if ζ is the rate at which such perturbation is transmitted between nodes, and if s i (t) and x i (t) are the probabilities that the node i is susceptible to be perturbed or get perturbed at time t, respectively, we can write the dynamics [161]:

dx i (t) dt = ζ (1 -x i (t)) j∈N A ij x j (t) , t ≥ t 0 , (83)
where A ij are the entries of the adjacency matrix of the graph for the pair of nodes i and j, and N is the set of nearest neighbors of j. In matrix-vector form the model becomes [161]:

ẋ (t) = dx (t) dt = ζ [I N -diag (x (t))] Ax (t) , (84)
with initial condition x (0) = x 0 . Let us call this model the Unperturbed-Perturbed (UP) one, in full analogy with the Susceptible-Infected (SI) model used in epidemiology.

The temporal evolution of the two states are illustrated in Fig. 26. Theorem 11 Consider the network UP (SI) model with ζ > 0. For a connected graph with adjacency matrix A, the following statements hold:

1. if x 0 ∈ [0, 1] n , then x (t) ∈ [0, 1] n , ∀t > 0.
Moreover, x (t) is monotonically nondecreasing i.e., x (t 1 ) ≤ x (t 2 ) for all t 1 ≤ t 2 . Finally, x 0 > 0n then x (t) ≫ 0n, ∀t > 0;

2. the UP model has two equilibrium points: 0n (no perturbation), and 1n (full perturbation); 3. the linearization of UP model about the equilibrium point 0n is ẋ (t) = ζAx (t) which is exponentially unstable; 4. each trajectory with initial condition x (0) ̸ = 0n converges asymptotically to 1n, that is, the perturbation spreads monotonically to the entire network.

The SP model can be rewritten as

1 1 -x i (t) dx i (t) dt = ζ j∈N A ij 1 -e -(-log(1-xj (t))) , (85)
which is equivalent to

dy i (t) dt = ζ j∈N A ij f (y j (t)) , (86)
where

y i (t) := g (x i (t)) = -log (1 -x i (t)) ∈ [0, ∞], f (y) := 1 -e -y = g -1 (y).
dŷ (t) dt = ζAdiag (1 -x (t 0 )) ŷ (t) + ζb (x (t 0 )) , (87)
where x (t) = f (ŷ (t)) in which x (t) is the approximate solution to the SP model, ŷ (t 0 ) = g (x (t 0 )) and b (x

) := x + (1 -x) log (1 -x) .
Let us then consider that a perturbation at a given node has a delay, a memory, before it is transmitted to any nearest neighbor. In this case we should replace the time derivative in the previous linearized model by a fractional time derivative, such that

[1] D α t ŷ (t) = ζ α Adiag (1 -x 0 ) ŷ (t) + ζ α Ab (x (0)) , (88)
where x (t) = f (ŷ (t)) in which x (t) is an approximate solution to the fractional SI model, ŷ is the solution of ( 88) with initial condition ŷ (0) = g (x (0)) and b (x) := x + (1 -x) log (1 -x) . Here D α t f (t) is the fractional time derivative in the Caputo formulation [START_REF] Caputo | Linear models of dissipation whose Q is almost frequency independentii[END_REF], namely

D α t f (t) := 1 Γ (⌈α⌉ -α) t 0 f (⌈α⌉) (τ ) dτ (t -τ ) α+1-⌈α⌉ . ( 89
)
For convenience, we write Ω := diag (1 -x 0 ) , and  = AΩ. It was then proved that this solution is an upper bound to the exact fractional SI model.

Theorem 12 [START_REF] Abadias | Fractional-order susceptible-infected model: denition and applications to the study of COVID-19 main protease[END_REF] For any t ≥ 0, we have

x(t) ⪯ x(t) = f (ŷ(t)) ⪯ x(t),
under the same initial conditions x 0 := x(0) = x(0) = x(0), where the solution ŷ of ( 88) is given by

ŷ (t) = E α,1 (ζt) α Â g (x 0 ) + ∞ n=0 (ζt) α(n+1) Ân Ab (x 0 ) Γ (α (n + 1) + 1) . ( 90
)
Furthermore, ∥x(t) -x(t)∥ → 0 and ∥x(t) -x(t)∥ → ∞ as t goes to innity.

Corollary 3 [START_REF] Abadias | Fractional-order susceptible-infected model: denition and applications to the study of COVID-19 main protease[END_REF] Let

x 0 ⪯ 1, x 0 = c N where c ∈ R + , let γ = 1 -x 0 . Noting that diag (1 -x (0)) = γI, then ŷ (t) = 1 -γ γ E α,1 t α ζ α γA ⃗ 1 - 1 -γ γ + log γ ⃗ 1, (91)
where E α,1 . is the Mittag-Leer matrix function of the corresponding matrix.

Therefore, here again if we take the sum of the entries of ŷ (t) at a given t we have

n i=1 ŷ (t) = C 1 tr E α,1 t α ζ α γA + tr (J -I) E α,1 t α ζ α γA -nC 2 , (92)
We recall that the best dierentiation obtained between the main proteases of Remark 10 The results found in this section indicate that the PRN of the main protease of SARS-CoV-2 has more capacity of transmitting perturbations across the whole protein than that of SARS-CoV-1, particularly to transmit such perturbations at longer distances in the protein.

Drug repurposing studies

The structure of the SARS-CoV-2 main protease bounded to an inhibitor has been determined crystallographically in several studies. For instance, Jin et al. [119] determined the structure of Mpro bounded to the inhibitor denoted as N3 and Zhang et al.

[260] determined the 3D structure of the main proteinase of SARS CoV-2 bounded to an α-ketoamide inhibitor. These structures have been used in several studies for investigating existing drugs that can be repurposed for targeting SARS-CoV-2 (see Table 10).

The main goal of these studies is to predict the best conformations/orientations of a ligand within the protein binding site. Docking process consists in eectively sampling the conformational space described by the free energy landscape to nd conformations/orientations that minimize a given scoring function, which should associate the native bound-conformation to the global minimum of the energy hypersurface. In Table 10 With the goal of improving the mechanisms of drug repurposing other studies have investigated other potential mechanisms and targets of SARS-CoV-2 main protease.

For instance Nguyen et al. [179] have identied 13 dierent binding sites in this protein.

However, the binding pocket, which is around the catalytic site is the most populated one by inhibitors. The only studies published so-far that use networks for analyzing the potential mechanism of inhibition of SARS-CoV-2 main protease were published by Estrada et al. [START_REF] Estrada | Topological analysis of SARS CoV-2 main protease[END_REF]. In these studies a few structures of the SARS CoV-2 protease bounded to three inhibitors, two of them being the most potent ones reported so far, were analyzed. The crystal structures of the three complexes have PDB codes: 6M0K

[52], 6LZE [START_REF] Dai | Structure-based design of antiviral drug candidates targeting the SARS-CoV-2 main protease[END_REF] and 6Y2G [260]. In these studies the author calculated the dierence:

∆G α ij = 1 n (n -1) i̸ =j G α ij (bounded) -G α ij (free) G α ij (free) , (93)
where G α ij = E α,1 (ζt) α γA ij . The parameters used in the study are: ζ = 0.01, γ = 1 -c n with c = 0.005. The results are given in Table 11. The values of the average path length L for paths between the top ten pairs of amino acids according to ∆G α ij and the number of times a residues in one of these paths is located in the binding site of the protease, N BS , are also given in the Table 11.

The results show better consistency for α = 0.5 where the most potent inhibitor, 6M0K, has the biggest eect in transmitting the perturbations from the binding site to the rest of the protein, i.e., the biggest ∆G α ij , but also accompanied by the largest α = 1.0 α = 0. This strategy developed by Estrada [START_REF] Estrada | Topological analysis of SARS CoV-2 main protease[END_REF] based on PRN and network descriptors was recently applied by Chen et al. [START_REF] Chen | Review of COVID-19 antibody therapies[END_REF] to the analysis of seven existing antibodies for SARS-CoV-2 spike (S) protein with three-dimensional (3D) structures deposited in the PDB. Five antibody structures associated with SARS-CoV were evaluated for their potential in neutralizing SARS-CoV-2.

Conclusion

We have shown how the combination of mathematical techniques provide physical meaning to network virology analysis and how these techniques are able to give important insights towards drug design and drug repurposing strategies.

Protein-protein interactions in virology

In general, viruses are small machineries which have very few proteins that interact with some of the proteins of their hosts [START_REF] Goodacre | Proteinprotein interactions of human viruses[END_REF]140,164]. This is resumed in Table 12 for some viruses of importance due to their impact on humans. None of these viruses has more than 100 proteins, with most of them having less than 20 proteins. This is a tiny fraction in comparison with the number of proteins estimated in humans, which ranges between 80,000 and 400,000 ones [247]. Another interesting characteristic of viral proteomes is the very small number of virus-virus (vv) protein-protein interactions (PPI). For instance, 50 proteins of the 89 one existing in the KSHV interact with each other through 123 PPI. Other viruses having more proteins, such as vaccinia virus which has 266 proteins, also have very few vv-PPI, i.e., 37 in the case of vaccinia. For the case of viruses with small number of proteins the number of vv-PPI is ridiculously small, e.g., 6 in HCV. Obviously, this number is already bounded by the small number of proteins as it cannot be larger than n (n -1) /2, where n is the number of proteins.

As can be seen in Fig. 27 once the virus enters a host cell and discharges its genetic material, its proteins can interact with those of the host, e.g., human, forming a series of virus-host (vh) PPI [START_REF] Goodacre | Proteinprotein interactions of human viruses[END_REF]140,164]. In Table 12 we report the number of viral proteins that interact with human ones (PIH) as well as the number of vh-PPI for the 7 viruses considered before. In average, every viral protein of the 7 viruses in Table 12 interacts with 24 human proteins. Such average is much bigger in the case of viruses like Zika (73.6) and Inuenzae A (84.1), while in others like KHSV (8.9) and HCV (13.9) are relatively small. Although, it is known that many PPI remain undetected for vh interactions the previous numbers reveal the magnitude of these important vh-PPI [START_REF] Goodacre | Proteinprotein interactions of human viruses[END_REF]140,164]. Table 12: Number of proteins of some important viruses as well as the number of them found to interact with human proteins (PIH) as well as the number of virus-human protein-protein interactions (VH-PPI).

Transmission of perturbations in virus-host PPI

The spread of perturbations across protein-protein interaction network is considered as one of the major causes of diseases [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF]206]. Such perturbations can be either of topological nature, e.g. deletion of nodes (proteins) or edges (interactions), or dynam-ical, i.e., the propagation of changes in the concentrations of given proteins in the cell.

One example of such perturbations is provided by the post-translational modication of host proteins by viral pathogens, which is used by the last to modulate host factors which are essential for viruses' replication, propagation, and evasion from host immune responses [197].

Nowadays it is possible to determine those human proteins which are targeted by viral proteins through PPI. We have seen in the previous section some examples of these VH-PPI networks. However, an important characteristic of these VH-PPI networks, which is frequently forgotten in the studies of the transmission of perturbations through these networks, is the following. The infecting viruses frequently invade the host cells at a specic organ, e.g., respiratory tract for respiratory viruses like respiratory syncytial virus, parainuenza viruses, metapneumovirus, rhinovirus, coronaviruses, adenoviruses, and bocaviruses. However, many of the proteins found in the VH-PPI are mainly expressed in organs outside the respiratory tract. In addition, there are evidence that support experimentally the importance of these extra-pulmonary VH-PPI. For instance, it is well-known that inuenza virus produces several important extrapulmonary damages [209], the most frequently ones being viral myocarditis and viral encephalitis, which indicates protein perturbations in the heart and brain of patients. In order to capture these important characteristics of the transmission of perturbations through a VH-PPI network we propose the mechanism and modeling strategy explained in the next subsection.

Mechanism

The mechanism described here was rst presented in [START_REF] Estrada | Protein-driven mechanism of multiorgan damage in COVID-19[END_REF][START_REF] Estrada | Cascading from SARS-CoV-2 to Parkinson's disease through protein-protein interactions[END_REF]. Here we continue its development by adding new details to the model nor presented before. We start by considering only those host proteins that are found to interact with viral proteins in the primary site of infection of the virus. That is, in the experimental settings for detecting VH-PPI, all viral proteins are interrogated to interact with all host proteins [START_REF] Goodacre | Proteinprotein interactions of human viruses[END_REF]140,164]. Therefore, some of the host proteins that are found to interact with viral ones are not necessarily expressed in the primary site of infection. These proteins can be important in cases of viral invasion of other organs where these proteins are mainly expressed. However, we consider here the initial stage of infection where the virus is mainly located at the primary site of infection.

The rst step in the transmission of perturbations is the determination of those host proteins in the primary site interacting with viral proteins which are able to perturb other host proteins independently on their location. For instance, let A be a host protein mainly expressed at the primary site which interact with a viral protein. Then, we ask if A is able to perturb another host protein B, which may be located at the primary site or not. We will clarify later what exactly we understand here by perturbing a protein.

If A does not perturb any host protein, it cannot be a propagator of perturbations and we simply discard it. However, if it perturbs B, then we ask for the location where B is mainly expressed. If B is mainly expressed in the primary site we consider that A can nd it through a diusive process and perturb it. Therefore, A is potentially a perturbator of B. However, if B is mainly expressed in an organ dierent from the primary site we need to ask whether the protein A can travel to this organ to trigger the perturbation of B. We will consider that this navigation from the primary site to the target organ is carried out as part of the cargo of an exosome (see later for explanation). Then, if A perturbs a protein B which is expressed outside the primary site, and A navigates outside the primary site as part of the cargo of an exosome, we consider that A is potentially a perturbator of B. This general mechanism is resumed in the owchart illustrated in Fig. 28.

Fig. 28: Flowchart developed here to identify potential protein perturbators of vulnerable proteins in humans.

Modeling strategy

1. A viral protein interacts with a host protein A mainly expressed in the primary site (PMEPS), which as a consequence, may be postranslationally modied (PTM)

[157, 252];

2. Let B be a protein also expressed in the primary site. Then, we consider that the PTM-PMEPS A navigates subdiusively in the primary site and may perturb B.

This process is modeled by considering the probability p B (t) that the protein B is perturbed at time t by a perturbation being propagated through the PPI network.

The evolution of p B (t) as the time progresses is determined by the fractional diusion equation (FDE) on the PPI subnetwork of PMELs (see Subsection 7.1.3):

D α t p B (t) = -C i∼B (p i (t) -p B (t)) , (94)
with a given initial condition (see further), where i ∼ B indicates that the sum is carried out over all nodes adjacent to B, and C is the diusion coecient of the system. The time-fractional derivative is dened for any 0 < α ≤ 1 and any function

f : [0, ∞) → R, as [156] D α t f (t) = 1 Γ (γ -α) t 0 f (γ) (τ ) dτ (t -τ ) α+1-γ ,
where Γ (•) is the Euler gamma function, γ = ⌈α⌉ and f (γ) represents the γth derivative of the function f ; 3. If the protein B is mainly expressed outside the primary site then we consider:

(a) whether the protein A is encapsulated in extracellular vesicles (EV), e.g., exosomes, and liberated from the primary site by exocytosis. This is an extremely complex set of processes involving numerous intracellular signaling reactions.

Therefore, it is considered here that if A has been reported to be part of the cargo of exosomes, then a xed excreted concentration of exosomes, f E P , is liberated by exocytosis at the primary site location;

(b) if A is encapsulated in exosomes we consider that it can be transported to the target organ (TO). This process is modeled by considering the time evolution of the number of exosomes containing a perturbed protein in its cargo (see Subsection 7.1.5);

(c) Once the PTM-PMEPS protein A is delivered at the TO, we consider that it navigates subdiusively and may perturb other TO proteins. This process is modeled again by using the FDE on the PPI subnetwork of the TO.

A graphical sketch of the general mechanism is illustrated in Fig. 29 where the nodes of a subnetwork of a hypothetical host PPI are split into two subsets, one representing proteins mainly expressed in the primary site, e.g., in the lungs, and the others representing proteins mainly expressed in an external organ, e.g., the brain. The connection between the two subnetworks is carried out by means of the transport of proteins in both directions mediated by exosomes.

We now explain the details of each of these steps in separate subsections.

Protein perturbations at the primary site

Once a PMEPS A is PTM by a viral protein it can navigate diusively inside the primary site and interact with other proteins also expressed in that site [START_REF] Batada | Stochastic model of proteinprotein interaction: Why signaling proteins need to be colocalized[END_REF]257]. It is well-known that the intracellular environment has a crowded nature with spatial barriers for the diusing protein [107], such that subdiusive processes more than normal diusion are expected for these protein-protein encounters. By subdiusive process we understand that the mean square displacement of a protein scales as [165] x 2 (t

) ∼ t κ , ( 95
)
where 0 < κ < 1 is the anomalous diusion exponent. Gupta et al. [107] used state-of-the-art neutron spin-echo (NSE) and small-angle neutron scattering (SANS), to study the diusion of two globular proteins in a crowded environment formed by poly(ethylene oxide) (PEO) which mimics a macromolecular environment. They found that in higher concentration of polymeric solutions, like in the intra-cellular space, the diusion is fractional in nature, which was mainly due to the heterogeneity of the polymer mesh in the bulk sample, which may well resemble the intra-cellular environment. Therefore, here we consider an FDE of the form [START_REF] Cao | Distributed coordination of networked fractional-order systems[END_REF]154,156]:

D α t x (t) = -CLx (t) , (96)
with initial condition x (0) = x 0 , where x i (t) is the probability that protein i is perturbed at time t, C is the diusion coecient of the network, which we will set hereafter to unity, and L is the graph Laplacian. Then, we have the following.

Theorem 13 The solution of the fractional time diusion model on the network is

x (t) = E α,1 -(tC) α L x 0 . ( 97
)
where E α,β (γL ) is the Mittag-Leer function of the Laplacian matrix of a graph.

Proof We use the spectral decomposition of the network Laplacian L = U ΛU -1 , where U = ⃗ ψ 1 • • • ⃗ ψn and Λ = diag (µr). Then we can write

D α t x (t) = -U ΛU -1 x (t) . (98)
Let us dene y (t) = U -1 x (t), such that D α t x (t) = -U Λy (t) and we have

U -1 D α t x (t) = -Λy (t) D α t y (t) = -Λy (t) . ( 99
)
As Λ is a diagonal matrix we can write

D α t y i (t) = -µ i y i (t) , (100)
which has the solution

y i (t) = E α,1 (-t α µ i ) y i (0) . ( 101
)
We can replace y i (t

) = U -1 x i (t) to have U -1 x i (t) = E α,1 (-t α µ i ) U -1 x i (0) x i (t) = U E α,1 (-t α µ i ) U -1 x i (0) , (102)
which nally gives the result in matrix-vector when written for all the nodes,

x (t) = E α,1 (-t α L) x 0 . (103)
We can write L = U ΛU -1 , where U = ⃗

ψ 1 • • • ⃗ ψn and Λ = diag (µr). Then, E α,1 (-t α L) = U E α,1 (-t α Λ) U -1 , (104)
which can be expanded as

E α,1 (-t α L) = ⃗ ψ 1 ⃗ ϕ T 1 E α,1 (-t α µ 1 ) + ⃗ ψ 2 ⃗ ϕ T 2 E α,1 (-t α µ 2 ) + • • • + ⃗ ψn ⃗ ϕ T n E α,1 (-t α µn) , (105)
where ⃗ ψ j and ⃗ ϕ j are the jth column of U and of U -1 , respectively. Because µ 1 = 0 and 0 < µ 2 ≤ • • • ≤ µn for a connected graph we have

lim t→∞ E α,1 (-t α L) = ⃗ ψ 1 ⃗ ϕ T 1 , (106)
where

⃗ ψ T 1 ⃗ ϕ 1 = 1. Let us take ⃗ ψ 1 = ⃗ 1, such that we have lim t→∞ ⃗ x (t) = lim t→∞ (E α,1 (-t α L)) ⃗ x 0 = ⃗ 1 ⃗ ϕ T 1 ⃗ x 0 = ⃗ ϕ T 1 ⃗ x 0 ⃗ 1.
(107)

Exosome-mediated transport of proteins

Exosomes are extra-cellular vesicles (EV), which are small, membrane-enclosed vesicles ranging from ∼ 40 nm to several microns in size, and which transport a cargo of proteins, lipids, metabolites and microRNA across organs [103,145,174,218,231].

EV have been found to play important roles in human diseases [START_REF] Akbar | Extracellular vesicles in metabolic disease[END_REF][START_REF] Alibhai | Cellular senescence contributes to age-dependent changes in circulating extracellular vesicle cargo and function[END_REF][START_REF] Frühbeis | Emerging roles of exosomes in neuronglia communication[END_REF]102,117,123,133,143,166,172,187,194,204,226,261,262]. More importantly, exosomes have been found to play fundamental roles in virus spread and replication [START_REF] Bello-Morales | Extracellular vesicles in viral spread and antiviral response[END_REF]137] particularly of HIV, HCV and SARS [START_REF] Giannessi | The role of extracellular vesicles as allies of HIV, HCV and SARS viruses[END_REF], Newcastle disease virus [264], HPV [104], avivirus infection [196], SARS-CoV-2 infection [110], as well as in the pathogenicity of certain viruses in the nervous system [138].

Here we consider a model that describes how the number of exosomes containing perturbed proteins evolves in time, such that we can estimate the approximate amount of these exosomes that arrives at a target organ. Let N E (t) be the number of exosomes containing a perturbed protein in its cargo. Then, the temporal evolution of N E (t) is

given by (see [225])

dN E (t) dt = E I (t) + E U (t) -k el N E (t) , (108)
where E I (t) and E U (t) are the rates at which exosomes produced by infected (I) and uninfected (U ) cells are released into the serum, and k el is the elimination rate of exosomes from serum.

The rate E I (t) is given by the product of the fraction f E I of exosomes containing perturbed proteins produced by infected cells entering the serum, the secretion rate R I of exosomes per infected cell (# /day/cell) and the total number N I of infected cells in the primary site that produce perturbed proteins at time t:

E I (t) = f E P R I N I . (109)
Let us consider that the rate E U (t) at which exosomes are produced by uninfected (U ) cells is negligible, such that we can assume E U (t) = 0.

Then, we consider that the number of infected cells in the primary site growth according to the Gompertz equation [225]:

N I (t) = N I (0) exp k growth k decay 1 -exp -k decay • t , ( 110
)
where N I (0) is the total number of infected cells in the primary site that produce perturbed proteins at time t = 0, k growth and k decay are, respectively, the growth rate to infected cell population (/day) and the rate at which abnormal growth rate decrease (/day).

Then, we have

dN E (t) dt = f E P × R I × N I (0) exp k growth k decay 1 -exp -k decay × t -k el × N E (t) . ( 111
)
We consider now some values for the dierent parameters in the model.

In the case of N I (0) we consider the number of infected cells in an individual with severe COVID-19. Then, using the values in Table 13 we calculate the evolution of the concentration of perturbed proteins in the blood: Table 13: Values of the parameters used in the model for exosome transport.

C (t) = N E (t)× NP (t) V T , ( 112
)
where NP (t) ≈ 10 3 is the average number of proteins per exosome [244] and V T ≈ 1.

Conclusions

1. the concentration of perturbed proteins may be very high even at very early times of infections, ranging between 10 7 -10 9 perturbed proteins in the blood of a patient with severe COVID-19 in the rst week. Therefore, here we will consider that such concentration is sucient to produce a signicant number of perturbed protein that penetrate a target organ and we will not simulate such process for individual proteins;

2. the peak of the concentration of perturbed proteins in the blood can be reached at later times of infection, e.g., more than 80 days, indicating a long persistence of these proteins in the blood of patients with severe COVID-19.

Extra-pulmonary eects of COVID-19

Many patients of COVID-19, particularly those with severe infection, display a series of multi-organ damages as a consequence of the disease [106,118,188,199,265]. Some of these damages reported in the medical literature are:

Cardiac: Potential perturbators: To identify the potential perturbators we proceed as follow. Let P be the set of all proteins in the human PPI network. First, we consider the set P I ⊂ P of all proteins which interact directly with SARS-CoV-2. We then identied those proteins in P I which are signicantly expressed in the lungs according to The Human Protein Atlas [235] and form the subset P L ⊂ P I . Then, we interrogate the database ExoCarta [126] for each protein in P L . Those proteins in P L found to be reported as part of any extracellular vesicle in ExoCarta forms the subset of proteins P E ⊂ P L ⊂ P I . Then, we pick a protein i ∈ P E and set p 0,i = 1 and p 0,j = 0 for all j ̸ = i where j is a node of the human PPI network. We then perform a simulation of the subdiusive propagation of a perturbation from the protein i ∈ P E to the rest using the FDE and its analytic solution. If for a given value of α and a given value of t < tc, the probability that a protein j ∈ V V is perturbed increases signicantly we consider that the protein i is a perturbator of j. We explored dierent values of 0 < α < 1 and observed that for small values of this parameter, e.g., for α ≤ 0.75, the increase in the probability of perturbation for some proteins occur instantaneously, which lack any physical sense. Then, we selected here α = 0.75 as the threshold and tc = 1 for all simulations. This process is repeated for each of the proteins in P E . This creates a new subset of proteins P C ⊆ P E which is formed by all those proteins targeted by SARS-CoV-2, which are signicantly expressed in the lungs, that can navigate outside this organ via EVs, and that can perturb other human proteins via subdiusive dynamics in the PPI network.

Vulnerable proteins: to identify vulnerable proteins we consider the set of proteins P P which are perturbed by any protein i ∈ P C . Then, for any protein j ∈ P P we interrogate version 7.0 of DisGeNET ([185]) which contains 1,134,942 genedisease associations (GDAs), between 21,671 genes and 30,170 diseases, disorders, traits, and clinical or abnormal human phenotypes, and 369,554 variant-disease associations (VDAs), between 194,515 variants and 14,155 diseases, traits, and phenotypes. We selected proteins with denitive or strong evidence of being involved in a human disease or syndrome. Then, we interrogated each of these candidates by using the disease network database created by Menche et al. [163],

which contains 299 diseases and 3173 associated genes. This database integrates disease-gene annotations from Online Mendelian Inheritance in Man (OMIM) and UniProtKB/Swiss-Prot with GWAS data from the Phenotype-Genotype Integrator database (PheGenI), using a genome-wide signicance cuto of p value ≤ 5×10 -8 .

The set of proteins P V ⊂ P P is formed by those proteins potentially involved in human diseases which are signicantly perturbed by pertubators, so we consider them the set of vulnerable proteins.

The results of this analysis is illustrated in Fig. 32 where it can be seen that we have identied 38 perturbators and a limited set of 26 vulnerable proteins (we explored a limited subset of the whole human PPI network) which are involved in 105 diseases. We then grouped these diseases according to the main organs/systems targeted by them.

This resulted in 13 categories: (1) neurologic, ( [START_REF] Althaus | Estimating the reproduction number of ebola virus (EBOV) during the 2014 outbreak in west africa[END_REF] metabolism, ( 12) respiratory, and ( 13) eye disorders [START_REF] Estrada | Protein-driven mechanism of multiorgan damage in COVID-19[END_REF].

Conclusions These results show that the proposed mechanism of transmission of perturbations through a PPI network is able to describe all the extra-pulmonary damages produced in patients with severe COVID-19. Also important is the fact that both, perturbators and vulnerable proteins, can be used as potential targets of drugs to palliate the extra-pulmonary damages of COVID-19 as illustrated in Fig. 33.

In order to nd such drugs we interrogate existing drug repository databases and found 27 drugs which are inhibitors of vulnerable proteins found [START_REF] Estrada | Protein-driven mechanism of multiorgan damage in COVID-19[END_REF]. These drugs can act over the vulnerable protein either indirectly, e.g., via the interaction with protein perturbators of such proteins (route I in Fig. 33) or directly on the vulnerable protein (route II in Fig. 33). The complete list of drugs found is illustrated in Fig. 34.

Viral network epidemiology

One of the main characteristics of the diseases produced by viruses is that they are transmissible from one individual to others, except for dead end host scenarios [START_REF] Flint | pathogenesis and control[END_REF].

Therefore, a snapshot of the disease at a given time cannot provide a global picture of the whole process due to the fact that the risk that an individual has of acquiring the infection changes dynamically as levels of infection change temporally in the population. This implies that we need transmission models for understanding the dynamics describing the temporal evolution of the spread of infectious agents in a given population. Such dynamical models would allow us to assess which public health interventions to apply to mitigate the disease and when it is the best time to apply them. One of the main goals of these dynamical models is to determine the number of infected individuals at a given time. This is important to determine the number of new infections arising per unit time, which is known as the incidence [START_REF] Flint | pathogenesis and control[END_REF]249]. It is usually expressed at x% per year, or x cases per 1000 per year, or x cases per 100 000 per year. This concept should be dierentiated from that of prevalence. Disease prevalence [START_REF] Flint | pathogenesis and control[END_REF]249] is the proportion of the population that is infected at a given time, i.e., the number of individuals infected at a given time divided by an appropriate measure of the population. The prevalence is usually expressed as a percentage. The contrast between both concepts is better illustrated with an example. A highly infectious and lethal disease has a high incidence because the number of new infections arising per day is Another important and necessary dierentiation is about the incubation period of a virus, the symptomatic and infectious periods as illustrated in Fig.

35 [START_REF] Flint | pathogenesis and control[END_REF]. The rst is the time between being infected and showing signs and symptoms of being ill. In this period of time most of the processes occur at the microscopic level as the virus is sequestering the host cellular machinery and the host immune system is amplifying its response to that invasion. The time from which the host has clear symptoms of infection until they clear is the symptomatic period. When the host is shedding enough infectious virus that other host are a risk of getting infected is called the infectious period. This is dierent from the incubation period. It could be a presymptomatic period in which the host are infectious such as in the case of measles in which the host is infectious some days before symptoms are evident. However, there are other viruses like Ebola in which the infectious period starts only after the symptoms appear.

For a viral epidemic to occur there are two necessary, although not sucient, ingredients. First, there should be part of the population which is susceptible to get infected by the virus. Second, there should be part of the population capable of infecting susceptible ones. Then, at early times the prevalence of infection rises as infection spreads. As a consequence there is an increase in the incidence of infection, which triggers an increase in the prevalence and the epidemic accelerates. However, at a given point the supply of susceptible individuals drops as they have become infected and the incidence falls, even when the prevalence may continue to rise for a time. This fall of incidence triggers a fall in prevalence due to the fact that infected individuals recover, die or move away from the studied population at a rate which is faster than that at which susceptible individuals become infected. The course of the infection depends on whether the virus is able to persist in the population for long periods of time, such as the disease become endemic or it will go extinct at this population.

We have said before that the existence of susceptible and infectious individuals was a necessary but not sucient condition for an epidemic to outbreak. The missing ingredient here is the mean number of new infections caused by a single infectious individual in the population of interest. This number is known as the eective reproduction number R (t) and it dictates whether there would be epidemic or not [236]. Basically, an epidemic requires R (t) > 1 in order for the prevalence of infection increases due to the fact that more than one new infection arises from the average infected host before numbers are specic to a particular infectious agent and susceptible population, but R (t) depends also on the particular time of the epidemic. When R 0 < 1 the number of infectious individuals decreases monotonically to 0 and when R 0 > 1 this number rst increases and then drops, such that R 0 = 1 is a sharp bound between a disease dying out or becoming an epidemic. Some typical values of R 0 for viral diseases are given in Table 14.

Disease outbreak and location R 0

Ref.

H1N1 inuenza in South Africa (2009) 1.33 [248] Ebola in Guinea (2014) 1.51 [START_REF] Althaus | Estimating the reproduction number of ebola virus (EBOV) during the 2014 outbreak in west africa[END_REF] Zika in South America (2015-16) 2.06 [START_REF] Gao | Prevention and control of Zika as a mosquito-borne and sexually transmitted disease: a mathematical modeling analysis[END_REF] SARS epidemic in (2002-03) 3.5 [105] Smallpox in Indian subcont. (1968-73) 4.5 [START_REF] Anderson | Infectious diseases of humans: dynamics and control[END_REF] Poliomyelitis in Europe (1955-60) 6 [START_REF] Anderson | Infectious diseases of humans: dynamics and control[END_REF] Measles in Ghana (1960-68) 14.5 [START_REF] Anderson | Infectious diseases of humans: dynamics and control[END_REF] Table 14: Some examples of the basic reproductive number of well-known epidemics suered between middle of XX century to nowadays.

Compartmental models

A usual way of modeling the transmission of a disease is to split the population into dierent compartments. Let us consider that we have a population which is divided into the class of susceptible individuals (S) and the infected one (I) as illustrated in the left-hand side of Fig. 36. Then, if we remove the imaginary barrier that separate both compartments we allow for a reaction to occur, such that susceptible individuals can be infected by infectious ones with a given probability. The rate at which susceptible individuals become infected is called the birth rate of infection and represented here by β. If we allow that the reaction is reversible and infected individuals can return to the susceptible state with a rate γ we are in presence of a model known as Susceptible-Infected-Susceptible (SIS) one:

Our main assumption here is that after removing the barrier separating the two compartments the population is well mixed, i.e., we have stirred the reactor to obtain homogeneity. In this scenario we have that the rate at which the population of infected individuals changes in time is given by:

İ (t) = βS (t) I (t) -γI (t) , (113)
which for a given initial condition I (0) = I 0 ∈ [0, 1] and for β > 0, γ > 0 and β ̸ = γ has the solution:

I (t) = (β -γ) I 0 βI 0 -e -(β-γ)t (γ -β (1 -I 0 )) . (114)
The main problem with this scalar (SIS) model is that it does not consider the real contact between individuals. That is, it is possible that as illustrated in Fig. 36 (bottom right) there are individuals with more contacts than others. From a network perspective the scalar model corresponds to the case where every individual has the same number of contacts, equal degree, which corresponds to a regular graph.

To capture the inuence of the network of connections in the propagation of the viral disease we describe the change of the probability of getting infected with time according to the following equation [161]

İi (t) = β (1 -I i (t)) n j=1 A ij I j (t) -γI i (t) , (115)
which in matrix-vector form is written as:

İ (t) = β (1n -diag (I (t))) AI (t) -γI (t) . (116)
We can now compare the eects of using a scalar model versus a network-based one.

We build a regular graph with n = 1000 nodes in which every node has 8 connections. In contrast we build a Barabási-Albert random graph with n = 1000 nodes and mean degree equal to 8, but in which the distribution of such degrees is skewed. Indeed, there is one node with degree larger than 200, while most of the nodes have degree one. We use β = 0.008 and γ = 0.01 and started the infection with only one node infected. In Fig. 37 we illustrate the dierences in the progression of the disease. As can be seen the consideration of the network structure implies that the epidemics growths very quickly reaching its peak at much earlier time than what is predicted by the scalar model. This is mainly a consequence of the existence of super-spreaders in the networked system, which are those nodes with a degree much larger than the mean degree of the population.

We now state one of the fundamental results concerning the SIS model [161].

Theorem 14 Let Eq. ( 116) describes an epidemics spreading over a simple connected (undirected) graph G. Let λ 1 be the spectral radius of the adjacency matrix A of G, and let β > 0, γ > 0. Then,

1. if I 0 , S 0 ∈ [0, 1] n , then I (t), S (t) ∈ [0, 1] n , ∀t > 0. Particularly, if I 0 > 0n then I (t) ≫ 0n, ∀t > 0;
2. 0n is an equilibrium point, and the linearization of Eq. ( 116) about it is (a) the linearization is exponentially stable (b) if I 0 ̸ = 0n then, t → ψ T 1 I (t) is monotonically and exponentially decreasing, and all the trajectories converge to 0n. 4. if βλ 1 /γ > 1, then (a) the linearization at 0n is exponentially unstable, i.e., there is an epidemic outbreak;

İ (t) = (βA -γ1n) I (t) ;
(b) there is a second equilibrium point I ⋆ , called the endemic state, such that i. If y (0) ≥ 0 is a scalar multiple of ψ 1 and satises either 0 < max i y i (0) ≤ 1 -γ/ (βλ 1 ) or min i y i (0) ≥ 1 -γ/ (βλ 1 ), then

I ⋆ ≫ 0n, ii. I ⋆ = (βλ 1 /γ -1) n i=1 ψ 3 1,i -1 ψ 1 +O (βλ 1 /γ -1) 2 as (βλ 1 /γ -1) → 0 + , iii. I ⋆ = 1n -(γ/β) K -1 1n + O γ 2 /β 2 as γ/β → 0 + ,
lim k→∞ y (k) = I ⋆ . (119)
Moreover, if max i y i (0) ≤ 1 -γ/ (βλ 1 ), then y (k) is monotonically non- decreasing, while if min i y i (0) ≥ 1 -γ/ (βλ 1 ), then y (k) is monotonically nonincreasing.

v. The endemic state I ⋆ is locally exponentially stable and its domain of attraction is [0, 1] n \0n.

Remark 11

The basic reproduction number for the network SIS model is R 0 = β λ 1 /γ.

This number indicates the threshold below which the disease dies out, while over it the disease becomes an epidemic.

Remark 12 When γ = 0 the infected individuals cannot return to the susceptible state.

In this case the SIS model is reduced to the simplest compartmental model, which is the SI one. This model was previously described for the case of the unfolded-folded transition of proteins.

Ṡi (t) = -βS i (t) n j=1 A ij I j (t) , İi (t) = βS i (t) n j=1
A ij I j (t) -γI i (t) , Ṙi (t) = γI i (t) .

(120)

In matrix-vector form it is written as

Ṡ (t) = -β (diag (S (t))) AI (t) , İ (t) = β (diag (S (t))) AI (t) -γI (t) . (121)
Theorem 15 [161] Let Eq. ( 121) describes an epidemics spreading over a simple connected (undirected) graph G. For t > 0 let ϱ 1 (t) and φ 1 (t) be the spectral radius and corresponding eigenvector of the matrix diag (S (t)) A , and let β > 0, γ > 0. Then, 1. if I 0 > 0n, and so S (t) ≫ 0n, then, (a) t -→ S (t) and t → I (t) are strictly positive ∀t > 0, (b) t -→ S (t) is monotonically decreasing, and (c) t -→ ϱ 1 (t) is monotonically decreasing, 2. the set of equilibrium points is the set of pairs S ⋆ , 0 , for any S ⋆ ∈ [0, 1] n . The linearization of the network SIR model about the equilibrium point S ⋆ , 0 is:

S (t) = -βdiag S ⋆ AI, I (t) = βdiag S ⋆ AI -βI; (122)
3. let βϱ 1 (τ ) < γ for τ ≥ 0, then t → φ 1 (τ ) T I (t) decays exponentially to zero for t ≥ τ ; 4. let βϱ 1 (τ ) > γ and I (t) > 0n, then (a) t → φ 1 (0) T I (t) grows exponentially fast with rate βϱ 1 (0) -γ, indicating an epidemic outbreak and (b) there exist τ > 0 such that βϱ 1 (0) < γ; 5. each trajectory converges asymptotically to an equilibrium point, that is, lim t →∞ I (t) = 0n, so that the epidemic asymptotically dies.

Remark 13

The eective reproductive number for the network SIR model is R (t) = β ϱ 1 (t) /γ. If we consider that very close to t = 0 the majority of individuals are susceptible, i.e., S (0) = 1, then R 0 = β λ 1 (A) /γ, where λ 1 (A) is the spectral radius of the adjacency matrix ([161]).

In several viral diseases there is a latent period in which an individual is infected but not yet infectious. In compartment models, this period is named Exposed (E) and it is an intermediate period between being susceptible (S) and becoming infectious (I). The resulting model is a four compartment one know as SEIR, whose scheme and equations for the scalar model are given below:

In the case in which a network of interactions is taken into account, the system can be expressed in matrix-vector form as [153]:

Ṡ (t) = -βdiag (S (t)) AI (t) , Ė (t) = βdiag (S (t)) AI (t) -σE (t) , İ (t) = σE (t) -γI (t) , Ṙ (t) = γI (t) .
R 0 = β β + γ E [k] + var (k) E [k] , (124)
where E [k] is the expected degree taken with respect to the degree distribution and var (k) is the variance of the degrees in the graph. Then, by keeping the same birth and death rate of an epidemic and keeping xed the mean degree of a network, the basic reproduction number increases linearly with the variance of the degree. In the scalar model var (k) = 0 meaning that R 0 is the minimum, and in networks with power-law degree distribution where var (k) is very large, the basic reproduction number is also large. The term β/ (β + γ) represents the probability that an infectious node transmits the infection to any of its nearest neighbors before it recovers from infection. A disease spreading becomes an epidemic if R 0 > 1, which means that E

[k] + var (k) E [k] > (β + γ) /β.
Therefore, for the same epidemiological conditions (right-hand side of the equation), the networks with the higher variance will have more chances of suering an epidemic outbreak than those with more homogeneous degree distributions.

The main dierences between the four models previously stated (SI, SIS, SIR and SEIR) are represented graphically in Fig. 38 where we have simulated the progression of an epidemic in an Erd®s-Rényi random graphs with n = 1000 nodes and m = 4000

edges and using the parameters stated in the Figure 38.

9 Network epidemiology in COVID-19

Due to the relevance and impact on the world population of the COVID-19 there have been an extremely large number of publications concerning epidemiological modeling of the disease. Thus, it is impossible to make an account of all of them here and we give a condensed resume in Table 15 of some of the most relevant ones from the point of view of the combined application of compartmental models and networks. [113] Some of the most important epidemiological parameters of COVID-19 have been studied in diverse works. They have been meta-analyzed, and we provide a resume of these results below. In Table 16 we report the values of the incubation time obtained from a series of meta-analysis of the data previously published in the scientic literature. In two of the meta-analysis the authors considered the pooled eect size calculated using xed eect and random eects models, one study reports only the random eect and two report the mean. The values of the incubation time range from 5.12 to 7.04 days, which in general point out to about one week of incubation of the virus. Another interesting variable also resumed in Table 16 is the serial interval. It is the duration between symptoms onset of an infector and that of one of its infectees in a transmission chain. The values could be negative if the infectee shows symptoms before than its infector. We provide two meta-analysis reporting both xed and random effects and two others reporting only the mean values. All in all, the serial period ranges from 4.98 to 6.49 days, showing again a short period of time between the appearance of symptoms in infectors and their infectees. Table 16: Resume of the results of several meta-analysis for the incubation period and serial period of COVID-19.

In another meta-analysis of published data Ahammed et al. [START_REF] Ahammed | Estimation of novel coronavirus (COVID-19) reproduction number and case fatality rate: A systematic review and meta-analysis[END_REF] studied the basic reproduction number R 0 obtained for COVID-19 in dierent regions of the world. In 14).

In the same Table we also provide the case fatality rate (CFR) of COVID-19 as meta-analyzed by Ahammed et al. [3]. CFR is the proportion of patients who die from COVID-19 compared with the total number of patients identied with the disease within a dened population of interest. As can be seen this parameter shows a large variability between dierent regions, with the smallest value for Oceania and the larger ones for Africa. Globally, the data reveals that about 3% of patients identied with COVID-19 die relative to those diagnosed with the disease. This CFR of COVID-19 is similar to that of > 2.5% reported for 1918 ue pandemic, which was been called the mother of all pandemics.

Finally, we would like to mention some of the eects produced by the non-pharmaceutical Here we consider the context of the transmission of an infection through the airports forming the air transportation network of a given country or even the global one.

Therefore, we consider a network G where the nodes represent the airports and two nodes are connected if there is a direct ight between them. if there is a ight from A to be B, it is always the case that there is also a ight from B to A. Thus, the network is undirected. The rst important question is which compartmental model to select to perform the analysis. In Fig. 40 we illustrate the typical situations which emerges in this scenario. For instance, as illustrated in Fig. 40 (a) an airport can have infected individuals (left panel) or not (right panel). In the rst case we say that the airport is infected (I ), while in the second case we say that the airport is susceptible (S). Now, if there is at least one ight from an infected to a susceptible airport we can have the following two situations. First, it is possible that some of the infected individuals from I travel to S. In this case, both airports become infected as illustrated in Fig. 40 (b). the second possibility is that all the infected passengers in I travels to S, in this case the second airport becomes infected while the rst recover the state of susceptible as illustrated in Fig. 40 To motivate the problem we analyze here a simple example. In Fig. 41 we illustrate a small graph that represents a hypothetical airport system, i.e., every node represents an airport and two nodes are connected if there is at least one ight between them. Suppose that at the initial stages of the propagation of an infection we are interested in canceling some connections between airports such that we mitigate the propagation of the disease through the whole system. Suppose that under the epidemiological conditions of this epidemic, e.g., a SIS model with β = 0.008 and γ = 0.01, and supposing that the epidemic can arrive at any airport with the same probability p = 1/n, we can obtain the time at which, let say 90%, of the whole network is infected. This time t ⋆ can be estimated by using the upper bound we have found previously, which is a worse case scenario. Then, our goal is to nd which edge(s) to cut to increase t ⋆ . Let us consider that we are interested in removing only one edge at the time, which implies canceling all ights between a pair of airports. Removal of the edge [START_REF] Abadias | Fractional-order susceptible-infected model: denition and applications to the study of COVID-19 main protease[END_REF][START_REF] Abergel | Giant viruses[END_REF] increases t ⋆ by 15% as well as removal of edge [START_REF] Moreno | Evaluation of the potential incidence of COVID-19 and eectiveness of containment measures in Spain: a data-driven approach[END_REF][START_REF] Alibhai | Cellular senescence contributes to age-dependent changes in circulating extracellular vesicle cargo and function[END_REF]. Therefore, they seem to be better candidates for the rst removal than edges [START_REF] Ahammed | Estimation of novel coronavirus (COVID-19) reproduction number and case fatality rate: A systematic review and meta-analysis[END_REF][START_REF] Akbar | Extracellular vesicles in metabolic disease[END_REF] or [START_REF] Abergel | Giant viruses[END_REF][START_REF] Ahammed | Estimation of novel coronavirus (COVID-19) reproduction number and case fatality rate: A systematic review and meta-analysis[END_REF] which increase by about 7%. But, as we will see below this is only part of the complex picture.

Before COVID-19 pandemic a signicant fraction of all global freight capacity was carried as cargo in the belly of passenger aircrafts. This kind of transport has increased signicantly during the pandemics. To have an idea of the capacity of this kind of transport let us mention that long-haul aircrafts like the Boeing 777-300EP and Airbus A350 can transport up to 16,000 kgs of cargo. Such cargo mainly consists of perishables, urgent goods, spare parts for land vehicles, spare parts for aerospace industry, drugs, vaccines and pharmaceutical products, machinery and accessories for medical use, among others. Therefore, removing the edge [START_REF] Abadias | Fractional-order susceptible-infected model: denition and applications to the study of COVID-19 main protease[END_REF][START_REF] Abergel | Giant viruses[END_REF] will make that transporting cargoes between these two airports will need 6 steps by using remaining connections. In fact, the average shortest ξ-path (as well as the average shortest topological distance)

increases by 33% in the network as a consequence of this single removal. This increase in the average distance needed to transport goods and passengers in the network represents a signicant increase in both time and cost for the whole system. In contrast, removing the edge [START_REF] Ahammed | Estimation of novel coronavirus (COVID-19) reproduction number and case fatality rate: A systematic review and meta-analysis[END_REF][START_REF] Akbar | Extracellular vesicles in metabolic disease[END_REF] that delays t ⋆ by 7%, does not change the shortest ξ-path distances and increases the shortest topological one by only 1%. Even the simultaneous removal of edges [START_REF] Ahammed | Estimation of novel coronavirus (COVID-19) reproduction number and case fatality rate: A systematic review and meta-analysis[END_REF][START_REF] Akbar | Extracellular vesicles in metabolic disease[END_REF] and [START_REF] Abergel | Giant viruses[END_REF][START_REF] Ahammed | Estimation of novel coronavirus (COVID-19) reproduction number and case fatality rate: A systematic review and meta-analysis[END_REF] increases only by 5% the shortest path distances with an increase of t ⋆ of 13%.

The previous simple example illustrates the complexity of the problem we are dealing with. Thus, we propose here to account simultaneously by the delay in the propagation of the epidemic and the relative increase in the average shortest ξ-path length.

Mathematical analysis

Let us start by considering the probability that the node v is not infected, which is

given by 1 -Iv (t). Then, the surprise or information content that the node v is not infected is given by [144]:

Sv (t) := -log (1 -Iv (t)) ∈ [0, ∞) , (125)
which makes that Ṡv (t) = İv (t) e Sv(t) . Therefore, this change of variable allows us to write the equation of the SIS model as follow [START_REF] Abadias | Fractional-order susceptible-infected model: denition and applications to the study of COVID-19 main protease[END_REF]: t) -γ e Sv(t) -1 .

Ṡv (t) = β w Avw 1 -e -Sv(
(126)

Now, using the facts that f (y) = 1 -e -y and g (y) = e y -1 are increasingly concave and convex functions, respectively, we have [START_REF] Bartesaghi | Where to cut to delay a pandemic with minimum disruption? mathematical analysis based on the SIS model[END_REF] Ṡv (t) < βq w Avw Ŝv (t)-

γ q w δvw Ŝv (t)+β (p + q log q) w Avw-γ p + log q q =: Ṡv (t) , (127)
where we have used I (0) = I 0 = p1n and let q = 1 -p, and where δvw is Kronecker delta function. In Matrix-vector form this equation is written as:

Ṡ (t) = βqA - γ q I Ŝ (t) + β (p + q log q) A - γ q (p + q log q) I 1n, (128)
which reveals that it is of the form:

Ṡ (t) = B Ŝ (t) + b. (129)
We now state a fundamental result about the solution of this approximate SIS model ( 128).

Theorem 16 [START_REF] Bartesaghi | Where to cut to delay a pandemic with minimum disruption? mathematical analysis based on the SIS model[END_REF] Let I (t), Ĩ (t) and Î (t) be, respectively, the solution of the exact, linearized, and approximate SIS model ( 128), with the same initial conditions: I (0) = Ĩ (0) = Î (0) = pIn. Then,

I (t) ⪯ Î (t) ⪯ Ĩ (t) , (130)
where a ⪯ b indicates that a i ≤ b i ,∀i ∈ V .

Let us now consider a bound for the solution of the approximate SIS model ( 128) for given initial conditions.

Lemma 7 [START_REF] Bartesaghi | Where to cut to delay a pandemic with minimum disruption? mathematical analysis based on the SIS model[END_REF] Let D := In -q 2 β γ A. Then, the surprise that a node v ∈ V in a graph G is not infected at a given time is bounded as

Iv (t) ≤ 1 -q exp   -C p q e qβtA 1n v -e γt/q e γt/q   , (131)
where

C =      max i x i if β γ > q 2 λ 1 -1 , min i x i if β γ < q 2 λ 1 -1 , (132)
with x i = I -pD -1 1n i .

Some of the general characteristics of the upper bound (131) for the SIS model are the following:

the term e qβtA 1n v contains all the structural information about the structure of the network. It represents the sum of all entries of the communicability matrix e ϱA corresponding to the node v, where ϱ = qβt;

the parameter C captures the epidemic threshold consistent with classic models; for small initial probabilities p it gives a close approximation of the exact solution and a very accurate description of the real spreading phenomenon, both above or under the threshold (see Fig. 42);

for large initial probabilities p, the solution still remains an upper bound, but the approximation gets worse as p → 1. Moreover, the bound predicts that Î (t) → 0 , below the threshold, only for small p.

Geometrization and shortest communicability paths

Let us assume that the SIS model is written in the form: Ṡ (t) = B Ŝ (t) + b, and let us simplify it to the form

Ṡ (t) = ζA Ŝ (t) , Ŝ (0) = I (0) , (133)
which has the solution Ŝ (t) = e ζA I (0) and where we remind that I (t) := 1 -e -ζS(t) . Let ev be the standard basis vector for node v, and let us make I (0) = ev. That is, we start an infection such that the probability that the node v is infected at time t = 0 is Iv (0) := 1 -e -1 and I w̸ =v (0) = 0. Then, we have that

Ŝv (t, ev) = e ζtA ev v , (134)
gives the probability that the node v is infected at time t, when the initial condition is I (0) = ev. So, as we have started with a very high probability that the node v is infected, the term Ŝv (t, ev) indicates the reinfection capacity of the node v.

we are accounting for the capacity of the node w of being infected by an infection which started at the node v. Similarly, Ŝv (t, ew) = e ζtA ew v accounts for capacity of the node v of being infected by an infection which started at the node w.

Denition 31 Let v and w be two nodes in a connected undirected graph. Then, the resistance of an infection to be transmitted between these two nodes is the sum of the reinfection capacities of each nodes, minus the sum of the capacities of each node of being infected from the other:

ξvw := Ŝv (t, ev) + Ŝw (t, ew) -Ŝw (t, ev) + Ŝv (t, ew) .

(136)

Remark 14 Notice that a low resistance to the infection transmission is possible if the two nodes have low capacities of retaining the infection by successive reinfections or if the capacities of each node of getting infected from the other is relatively high. Pairs of nodes with low resistance to transmission of an infection are the ones at the higher risk in a network.

We now have the following.

Theorem 17 Let v and w be two nodes in a connected undirected graph. Then, the resistance of an infection to be transmitted between these two nodes is given by: ξvw = (ev -ew) T e ζtA (ev -ew) ,

and it is a Euclidean distance between the two nodes.

Proof It is straightforward to realize that ξvw = (ev -ew) T Ŝ (t, ev) -Ŝ (t, ew) ,

Then, by using e ζtA = U e ζtΛ U T and having that φv is the vth column of U T , the last expression can be rewritten as

ξvw = e Λ/2 φv -e Λ/2 φw T e Λ/2 φv -e Λ/2 φw = (xv -xw) T (xv -xw) = ∥xv -xw∥ 2 , (140)
where xv := e Λ/2 φv is the vector of Euclidean coordinates of the vertex v.

The quantity ξvw gives a metric between every pair of nodes in a network. It measures how easy is the transmission of an infection from one node to another. However, it does identies directly the path followed by the infection from one node to another.

To understand this more clearly let us consider the illustration in Fig. Technically what we need is to make a geometrization of the graph such that we can identify the paths with the shortest sum of the distances ξvw along the edges of the path [START_REF] Bridson | Metric spaces of non-positive curvature[END_REF]158]. This is exactly what we do in the following.

Denition 32 Let us consider e = (v, w) as a compact 1-dimensional manifold with boundary ∂e = v ∪ w. Let the edge evw = (v, w) be given metric Lvw such that evw ∼ = isom [0, ξvw (ζt)] (v, w) ∈ E, 0 (v, w) / ∈ E.

(141)

We now extend the metric Lvw on the edges of G via inma of lengths of curves in the geometrization of G. Then, the network becomes a metrically length space, which is locally compact, complete and geodetic. We then have the following.

Denition 33 Let P (v, w) be a path in G connecting the nodes v, w ∈ V . Then, the shortest ξ-path length between the nodes v, w ∈ V is Cvw (G, ζt) := min

P (vw) (i,j)=e∈E e∈P ξ ij (ζt) , (142)
where the minimum is taken among all paths connecting the corresponding pairs of vertices.

Example 3 In the Fig. 44 we illustrate the shortest ξ-paths between the same pair of nodes in a random geometric graph with 400 nodes. As can be seen when ζt is relatively large the ξ-path of minimal resistance to the transmission of an infection between the two nodes goes by a set of nodes of relatively low degree. As soon as ζt drops the shortest ξ-path approaches to the shortest topological path between the two nodes.

Indeed, we have previously proved the following [216].

Lemma 8 Let ξvw (γ) = (exp (γA)) vv + (exp (γA)) ww -2 (exp (γA)) vw for a given γ > 0. Let Pvw (γ) be the shortest ξ-path between the nodes v and w for for a given γ > 0 and let Svw be the shortest (topological) path between the same pair of nodes.

Then, lim γ→0

Pvw (γ) = Svw.

(143)

Canceling ights with minimal disruption

Here we propose an algorithm for performing edge-removal in a real-time for a network in which an epidemic is propagating according to a SIS mechanism, such as in the case of an epidemics propagating through airports. We rst start by the following [START_REF] Bartesaghi | Where to cut to delay a pandemic with minimum disruption? mathematical analysis based on the SIS model[END_REF].

where the maximum is obtained among all the edges of the graph.

At a given value of t the index O (e, t) identies the edge e whose removal delays the most the epidemic propagation through the whole network with the minimal impact on the network navigability. Then, we design the following method of systematic edgeremoval on the network:

1. At time t 1 identify the edge e 1 that maximize O (e 1 , t 1 );

2. Remove the edge e 1 from G and model the evolution of the SIS epidemics for t > t 1 without the edge e 1 ; 3. At a given time t 2 repeat 1 and 2 until a given criterion of stop is reached. We apply the edge-removal strategy described in this work for removals at t = 200, 400, 600, 800, 1000. We normalized the edge weights by dividing the weighted adjacency matrix by the mean value of the edge weights in the network. We consider the edge removal strategy described before by considering the time at which 90% of the airports were infected. This percentage can be changed and dierent results can be obtained accordingly, but such exhaustive study is out of the scope of this book chapter.

The results of the edge removal strategy are illustrated in Table 18. By removing all the ights between Heathrow and Edinburgh the time for infecting 90% of the whole network drops by 38.9%. The second removal increases this time to 88.9% and the third one increases it up to 194.4%. Finally, the cumulative removal of 4 connections increases the time to infect 90% of the whole network by 333.3%. At the same time the remove of all ights between Heathrow and Edinburgh does not drop the capacity of the global network to diuse items/goods and passengers through its nodes. In contrast, it increases this capacity by 9.9%. This is, of course, a consequence of considering that such items/goods and passengers move in the network in a completely diusive way.

If we consider that they move using the shortest paths, then we observe a drop in the capacity of the network of only 0.05%, i.e., the increase in the average shortest path length in the network after the removal. After the four removals previously described the network has increased its diusive capacity by 3.6% with a drop in its capacity to deliver items/goods/passengers via shortest paths of 0.4%. In either way, the removal of these four inter-airport connections produces a remarkable delay on the propagation of the SIS disease in comparison with a very small aection of the network operative capacity [START_REF] Bartesaghi | Where to cut to delay a pandemic with minimum disruption? mathematical analysis based on the SIS model[END_REF]. It is important to remark that we have make the removals one at a time, such that the course of the epidemic is changed at every time an edge is removed. This is illustrated in Fig. 46, where it can be seen how the time evolution of the epidemics is changing continuously by modifying the topology of the network. Therefore, this emulates a real-world situation in which the progression of the epidemics is followed in real-time and such cancellations of origin-destination edges are done accordingly. Although this is intended here more like an academic exercise than as a real-world scenario, it is clear that the implementation of this strategy into more complex epidemiological models could be appropriate for delaying epidemic progression with minimal interruption of the economy.

11 Conclusions. The future of networks virology Network models are mathematical models. As in every model they make assumptions over which the description of a given phenomenon is approached. The goal of these models is to understand the mechanisms governing these systems and the implementation of potential scenarios that can be simulated on a computer. But, as every model, they are constrained by what it has being assumed on the rst place. When referring to models related to viral infections many authors focus only their attention on epidemiological models. That is, those models designed to forecast the evolution of an epidemic, as the compartmental models described here. James et al. [116] have stated that as a consequence of the proliferation of models, which frequently diverge in their conclusions, there have been criticism on the validity of modeled analyses and uncertainty as to when and to what extent results can be trusted. They remarked that this limitation has been mainly a consequence of basing such forecasts on a single model, either by using xed parameter values or by averaging the results of multiple parameter sets. In a Perspective paper in New Eng. J. Med., Holmdahl and Buckee [112] has pointed out that the accuracy of these models is constrained by our knowledge of the virus. In the case of an emerging disease like Covid-19, many biologic features of transmission are hard to measure and remain unknown.

Epidemiological models have been beneted from the advances in network theory. It is nowadays demonstrated beyond any reasonably doubt that considering the network of interactions is vital for a better predictability of an epidemic progression. However, still a large proportion of the models used for forecasting epidemics, like in COVID-19, are based on the assumption that the population is will-mixed and use scalar models (see [START_REF] Estrada | COVID-19 and SARS-CoV-2. modeling the present, looking at the future[END_REF] and references therein). It is true that knowing such contact networks in the real-world is a very dicult task. Then, it is necessary to investigate on how to infer these networks in a more realistic scenario than the classical models of network construction that exists today. Additionally, as we have seen here some of the biologic features needed to fed the epidemiological models can also be modeled and estimated by using mathematical network theory. Then, an obvious future area of development in networks virology is the integration of the dierent scales at which we can obtain information that can feed other levels of the whole process.

Some of the particular areas of networks virology are more developed than others.

In the area of protein residue networks, more interaction and integration with existing bioinformatic tools are necessary for advancing this subarea. Although many advances have been made on the area of protein-protein interaction networks, it remains unexplored how eects are transmitted between proteins which are in the same network but located at dierent, frequently very distant, organs. Here we have advanced a hypothesis which works well under its assumptions for explaining the extra-pulmonary eects of COVID-19. More works on this area in other dierent virus infecting humans are needed in combination with experiments that falsify this hypothesis.

We just will nish with a quotation made by George A. Bartholomew [23]: The complexity of contemporary biology has led to an extreme specialization, which has inevitably been followed by a breakdown in communication between disciplines. Partly as a result of this, the members of each specialty tend to feel that their own work is fundamental and that the work of other groups, although sometimes technically ingenious, is trivial or at best only peripheral to an understanding of truly basic problems and issues. There is a familiar resolution to this problem but it is sometimes diculty to accept emotionally. This is the idea that there are a number of levels of biological integration and that each level oers problems and insights that are unique to it;

further, that each level nds its explanations of mechanism in the levels below, and its signicance in the levels above it. Therefore, networks virology is a promising new area of research where much is still to be done, and where the future should point towards an interdisciplinary integration of contemporary virology and epidemiology in the context of complex systems analysis.

Fig. 1 :

 1 Fig. 1: Some of the major epidemics produced by viruses since 1880, according to data published by Piret and Boivin.

Fig. 2 :

 2 Fig. 2: Main geometrical shapes of viruses.

  : I: double stranded DNA viruses, e.g., Adenoviruses, Herpesviruses, Poxyviruses, Polyomavirus (SV40), Papillomavirus, Cytomegalovirus (CMV), Pox viruses (vaccinia) T4, T7 bacteriophages II: single stranded DNA (+) sense, e.g., Parvoviruses; III: double stranded RNA, e.g., Reoviruses; IV: single-stranded RNA viruses of positive (+) sense, e.g., Picornaviruses, Coronavirus, Hepatitis A/C/E viruses, Poliovirus, Rubella virus, Tobacco mosaic virus V: single stranded RNA viruses of negative (=) sense, e.g., Ebola, Marburg, Lyssavirus (Rabies), Morbillivirus (Measles), Rubulavirus (Mumps), Inuenza virus A/B/C, Lymphocytic choriomeningitis (LCMV) VI: positive (+) sense single stranded RNA viruses that replicate through a DNA intermediate, e.g., Retroviruses; VII: double stranded DNA viruses that replicate though an single stranded RNA intermediate, e.g., Hepatitis B virus. 2.1 SARS-CoV-2 Since December 2019 the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) has produced an outbreak of pulmonary disease which has soon became a global pandemic, known as COronaVIrus Disease-19 (COVID-19).

Fig. 4 :

 4 Fig. 4: Genomic classication of viruses according to Baltimore.

Fig. 5 :

 5 Fig. 5: Illustration of the life cycle of the SARS CoV-2.

Fig. 8 :

 8 Fig. 8: Examples of networks found in networks viroloy: (a) the protein residue network of SARS-CoV-2 spike S protein, (b) the protein-protein interaction network of SARS-CoV-2 (red circles) and human (blue circles) proteins, (c) a population contact network is the city of Haslemere, Surrey, U.K. (d) the network of airline connections between airports in the USA.

Fig. 9 :

 9 Fig. 9: Illustration of two isomorphic networks representing the amino acid interactions in the proteins with PDB code 1PPT.

Fig. 10 :

 10 Fig. 10: Illustration of the change in the dihedral angle φ between four consecutive atoms in a linear chain. The rst angle corresponds to a fully-folded structure and the last one to a fully-unfolded one.

Fig. 12 :

 12 Fig. 12: Part of a PDB le which contains the structural information of every atom in a protein.

  2H2Z [254], 2DUC [171], 1UJ1 [255], and 2BX4 [227]. For SARS CoV-2 we select the structures with PDB codes: 6M03 [258], 6M2Q [258], and 6Y2E [259] (see Table

Fig. 13 :

 13 Fig. 13: Illustration of the alignment of the amino acid sequences of 2BX4 and 6Y2E.

Fig. 14 :

 14 Fig. 14: Superposition of the 3D structures of the main proteases of SARS-CoV-1 (2BX4) and of SARS-CoV2 (6Y2E).

Fig. 15 :

 15 Fig. 15: Illustration of a linear chain of seven segments where the dihedral rotation is allowed (a). In (b) and (c) we have folded the planes around the edges of the segment to embed the structure into a cube. (d) and (e) The resulting linear chains folded inside a cube.

Fig. 16 :

 16 Fig. 16: Illustration of the potential function used for dening the folding degree index.

  (a) for the two structures determined at the lowest and the highest temperatures, respectively. Therefore, we should expect a drop of the I 3 index as the temperature rise as observed in Fig.17

  (b) where the Pearson correlation coecient is -0.92.

Fig. 17 :

 17 Fig. 17: (a) Superposition of the structures of the ribonuclease-A determined at 98 (red) and at 320K (blue). (b) Linear correlation between the temperature at which the structure of ribonuclease-A was determined and the corresponding normalized I 3 index.

Fig. 18 :

 18 Fig. 18: Illustration of six dierent conformations of a linear chain.

Fig. 19 :

 19 Fig. 19: Plot of the amino acid contributions to the folding degree for four crystallographic versions of the main protease of SARS-CoV-1 (a) and three versions of the one of SARS-CoV-2 (b).

Fig. 20 :

 20 Fig. 20: Dendrogram of the similarities between the amino acid contributions to the folding degree of the 4 structures of the main protease of SARS-CoV-1 and the 3 ones

Fig. 21 :

 21 Fig. 21: Illustration of the α-carbon backbone of the small protein with PDB code 1ppt (left), a visualization of the pairs of atoms at no more than 7.0 Å (center) and the resulting protein-residue network (right).

Fig. 22 :

 22 Fig. 22: Illustration of the small subgraphs appearing in Theorem 3.

.

Fig. 23 :

 23 Fig. 23: Dendrogram for the seven main protease structures considered in this Chapter based on the subgraph ngerprints analyzed here.

Fig. 24 :

 24 Fig. 24: Bar plot of the relative dierence (in percentage) between the average number of subgraphs for the three PRNs of SARS-CoV-2 protease minus those four of SARS-CoV-1. Notice that SARS-CoV-2 protease has about 10% more F 21 than SARS-CoV-1.

Remark 8

 8 Two examples of the use of this bound are (i) considering a = -λ 1 and b = -λn; (b) considering a = -kmax and b = kmax.

Table 7 :

 7 Mean and standard deviation of the values of the descriptors studied here for the four PRNs representing the main protease of SARS-CoV-1 and the three ones representing the main protease of SARS-CoV-2. The relative dierences between the means for the structures of both proteases are given as well as the probability that both samples comes from populations with the dierent variance according to the Ftest and the Kolmogorov-Smirnov test, respectively. Signicant values of the probability are bolded for p > 0.95.

6. 2 . 5

 25 How to account for longer walks? We now return to the problem of how longer subgraphs may be relevant in the differentiation of the main proteases of SARS-CoV-1 and SARS-CoV-2. Let us start by considering the four bicycles illustrated in Fig. 25. They correspond to the combination of chordless cycles of length 3 and 9 (I), 4 and 8 (II), 5 and 7 (III) and two of length 6 (IV).

Fig. 25 :

 25 Fig. 25: Illustration of four networks with the same number of vertices and edges but having dierent chordless cycles (holes): (I) a triangle and a nonagon, (II) a square and a octagon, (III) a pentagon and an heptagon, (IV) two hexagons.

α

  for the four crystallographic versions of the main protease of SARS-CoV-1 (top four) and the three versions of the SARS-CoV-2 one (bottom 3 structures).

Fig. 26 :

 26 Fig. 26: Illustration of the evolution of the reaction transforming unperturbed (U) (green curve) amino acids in perturbed (P) (magenta curve) ones according to the UP-model described in the text.

  Lee et al. [144] have considered the following linearized version of the previous nonlinear equation

  SARS-CoV-1 and SARS-CoV-2 was obtained by using indices based on Mittag-Leer functions of the adjacency matrix. Therefore, the Mittag-Leer (ML) subgraph centrality and communicability indices indicate the capacity of a given node in a PRN to transmit a perturbation across the whole network, such that it returns to the initial node (subgraph centrality) or ends up in a dierent one (communicability). So we have the following.

α

  for the three complexes of the SARS-CoV-2 main protease and inhibitors. The values of the inhibitory potency of each of the three inhibitors is given in the last column.values of L and N BS . On the contrary, the least potent inhibitor, 6Y2G, shows a decrease in the sensibility of the main protease to transmit perturbations across the structure, accompanied by the lowest values of L and N BS .

Fig. 27 :

 27 Fig. 27: Schematic representation of the virus-host protein-protein interactions generated after the entrance of the virus in the host cells.

Fig. 29 :

 29 Fig. 29: Scheme illustrating the transmission of perturbations across a PPI network of host where proteins are mainly expressed at dierent organs of the host. Such interorgan transport is illustrated here by means of extracellular vesicles, such as exosomes.

Fig. 31 :

 31 Fig. 31: Evolution of the number of exosomes liberated by infected cells at the primary site of infection assuming 10 4 (a) and 10 6 (b) infected cells. The rest of parameters are the ones provided in Table 13, NP (t) = 10 3 [244] and V T = 1.5 L [245].

Fig. 32 :

 32 Fig. 32: Interorgan perturbations through protein-protein interactions. Left column are proteins mainly expressed in the lungs with high probability of perturbing those with strong evidence of being involved in human diseases (central column), which are mainly expressed in the tissues given in the right column. The thickness of the connections is proportional to the number of diseases of a given class in which a protein is directly (2nd and 3rd columns) or indirectly (1st and 2nd columns) involved.

Fig. 33 :

 33 Fig. 33: General mechanism of inter-organ cascade giving rise to protein perturbations in COVID-19 patients which may be one of the causes of multi-organ damages in this disease. There are two points for potential drug interventions, which are marked as (I) for interrupting protein perturbators and (II) for interrupting interactions with vulnerable proteins.

Fig. 34 :

 34 Fig. 34: Repurposed drugs for interorgan damage. Group of known drugs (left column) which inhibit receptors mainly expressed in the lungs (second column from left) which have been found here to perturb proteins involved in human diseases (third column from the left) which are mainly expressed in several tissues or systems outside the lungs (right column).

Fig. 35 :

 35 Fig. 35: Illustration of the dierences between incubation period and symptomatic period during the course of an infection. From the epidemiological point of view the infectious period is the one in which an individual may propagate the virus to others and may coincide with the symptomatic period or extend to part of the incubation one.

Fig. 36 :

 36 Fig. 36: Illustration of the compartmentalization of a population into two classes: Susceptible (S) and Infected (I) (left panel). From the point of view of modeling the disease spreading (right panels) it can be considered that the population is well-mixed and use scalar models (top panel) or that a network of contacts exists (bottom panel). The results are signicantly dierent in both cases.

Fig. 37 :

 37 Fig. 37: Illustration of the evolution of the number of infected individuals in a population of 1000 individuals considering well-mixed population (regular) or a power-law distribution of the number of contacts per individual.

  where K is the diagonal matrix of node degrees. iv. Let us dene a sequence {y (k)} k∈N ⊂ R n by y (k + 1) := F + β γ Ay (k) .

A

  dierent model divides the population into three compartments: Susceptible (S), Infectious (I) and Recovered (R) which are those individuals that have recovered from the disease creating immunity to it or that have died. Then the Susceptible-Infected-Recovered (SIR) model has the following scheme and scalar equations [161]: In a network of interactions the SIR equations are [161]:

  importance of the skewness of the degree distribution of a network on the evolution of a SEIR model let us consider the following. First we consider the conguration model for generating random networks from a given degree sequence. In this context the basic reproductive number of the SEIR model is[153]

Fig. 38 :

 38 Fig.38: Evolution of the epidemic spreading using SI (a), SIS (b), SIR (c) and SEIR (d) models on an Erd®s-Rényi random graph with n = 1000 nodes and m = 4000 edges using β = 0.008, γ = 0.01 and σ = 0.02 for the respective models.

  because Ŝ (t, ev) = e ζtA ev and Ŝ (t, ew) = e ζtA ew by substituting in the previous equation we have ξvw = (ev -ew) T e ζtA (ev-ew)

43 .

 43 As can be seen the distance ξvw between the two nodes is equivalent to what an eagle (a) has to ight to go from one top of a hill to the other, but the mountaineer (b) has to go through the geodesic described by the paths going down the rst hill and then up the second. In a network it is impossible to ight like an eagle and it is necessary to nd the shortest path through the nodes and edges that goes from one node to another, like the mountaineer.

Denition 34 Let

 34 C (G, ϱ) be the mean shortest communicability path (SCP), that is the average of Cvw (G, ϱ) over all shortest diusive paths connecting pairs of nodes in G. Let G -e be the graph G from which the edge e has been removed. Then, the mean capacity of a network to reroute goods/items/passengers after the removal of an edge e is ∆ C (G -e, ϱ) := C (G -e, ϱ) -C (G, ϱ) t ⋆ be the minimum time at which a given percentage of the nodes in G are infected in a SIS dynamics on the graph. Let∆ t ⋆ (G -e, ϱ) := t ⋆ (G -e, ϱ) -t ⋆ (G, ϱ) t ⋆ (G, ϱ) . dene the following ratio as an indicator of the best candidate edge(s) to be removed,O (e, t) := max e∈E ∆t ⋆ (G -e, ϱ) ∆ C (G -e, ϱ) ,

10. 2 . 1

 21 UK airports networkIn Fig.45we show a representation of the network of domestic ights between 44 commercial airports in the United Kingdom in the year 2003. The network consists of 236 weighted edges representing air internal routes between these 44 airports. The weights correspond to the number of passengers transported during that year between the corresponding airports. The nodes are represented with size and colors proportional to their weighed degreestotal number of passengers arriving/departing to/from that airport in 2003. We remind the reader that 2003 was the year in which the SARS epidemic was spreading across the world.

Fig. 40 :

 40 Fig. 40: (a) Illustration of two airports with their respective populations. Because in the population of the rst airport there are infected individuals, this airport if considered to be infected (I) while the second is susceptible as there are no infected individuals in it. (b) As a result of a ight connection between the two airports it is possible that the population of infected individuals is split between the two airports. This results in the two airports becoming infected (I). (c) As a result of a igh connection between the two airports it is possible that the whole population of infected individuals moves tothe second airport, which will make that it becomes infected. In its defect the originally infected airport has returned to the susceptible (S) state, which is why the SIS model is the most appropriate for capturing this situation.

Fig. 41 :

 41 Fig. 41: Illustration of a toy graph used as an example of the edge removal process described in the text.

Fig. 42 :

 42 Fig.42: Evolution of the proportion of infected nodes in a SIS dynamics on a small graph with 12 nodes using the exact solutionIv (t) (black continuous curve), linearized solution Ĩv (t) (blue broken line) and the upper bound (131) Îv (t) (red dotted line). The parameters used are β = 0.004, γ = 0.001, p = 1/12. The largest eigenvalue of the adjacency matrix of this graph is λ 1 ≈ 3.098.

Fig. 43 :Fig. 44 :

 4344 Fig. 43: Illustration of the dierences between the ξ-distance between two nodes in a network (a) and the length of the shortest ξ-path connecting them (b). Because noninfectious particle can ight directly from v to w as can be done by an eagle in (a)it is necessary a geometrization of the graph to nd Lvw which would the route most probably followed by this infectious particle between the corresponding two nodes.

Fig. 45 :

 45 Fig. 45: Illustration of the UK airport connections network in 2003. The nodes have radii and are colored proportionally to the number of passengers that visited that airport in 2003.

Fig. 46 :

 46 Fig. 46: Change in the temporal evolution of an epidemic spreading through the airport connection network of the UK as a result of the edge-removal exercise implemented here.

Fig. 47 :

 47 Fig. 47: Schematic representation of the integration of dierent areas of network virology such that the output of certain areas feeds the parameters of other areas.

  The complement of a graph G is the graph Ḡ with the same set of nodes as G but two nodes in Ḡ are connected if and only if they are not connected in G. An empty or trivial graph is a graph with no links. It is denoted as Kn as it is the complement of the complete graph.

1. A forest is a disconnected graph in which every connected component is a tree.

A k-regular graph is a graph in which all nodes have degree k. A particular case of regular graph is the complete graph Kn previously dened, where every node has degree n -1. Another type of regular graph is the cycle, which is a regular graph of degree 2, i.e., a 2-regular graph, denoted by Cn.

Table 2 :

 2 List of the 20 naturally occurring amino acids, and their corresponding 3-and 1-letter codes. Denition 12 Let A 1 , A 2 , A 3 , A 4 be four consecutive atoms in the backbone of a protein P. The dihedral angle φ is dened as the angle formed between the planes

	Amino acid	3-letter	1-letter	Amino acid	3-letter	1-letter
	Alanine	Ala	A	Leucine	Leu	L
	Arginine	Arg	R	Lysine	Lys	K
	Asparagine	Asn	N	Methyonine	Met	M
	Aspartic acid	Asp	D	Phenylalanine	Phe	F
	Cysteine	Cys	C	Proline	Pro	P
	Glutamic acid	Glu	E	Serine	Ser	S
	Glutamine	Gln	Q	Threonine	Thr	T
	Glycine	Gly	G	Tryptophan	Trp	W
	Histidine	His	H	Tyrosine	Tyr	Y
	Isoleucine	Ile	I	Valine	Val	V

Table 3 :

 3 Proteins in the SARS CoV-2, their names and brief description of their func-

	tion.
		SARS CoV-1			SARS CoV-2
	PDB	res. (Å)	length	PDB	res. (Å)	length
	2H2Z	1.60	306	6M2Q	1.70	305
	2DUC	1.70	306	6Y2E	1.75	306
	1UJ1	1.90	301	6M03	2.00	306
	2BX4	2.79	298
	Table 4: Protein Data Bank (PDB) codes for structure of the main protease of SARS
	CoV-1 and SARS CoV-2 without inhibitors (apo forms).

  Denition 17 Let G = (V, E) be a graph in which the set of vertices represent the Cα atoms of a protein. Let two nodes of G be connected if and only if the inter-residue distance is equal or less than r C . This network corresponds to a PRN of the given protein. The adjacency matrix A of the PRN is then built with elements dened by

Table 8 :

 8 Values of Tr (E α,1 ) and Ῡ for two dierent values of the fractional parameter

Table 9 :

 9 5 both statistical tests, F -test and KS-test, dierentiate both samples as coming from populations with dierent variances. Indeed, the subgraph centrality and communicability for the PRN of SARS-CoV-2 are about 2,000% higher than those of SARS-CoV-1, which indicates that the main protease of SARS-CoV-2 is signicantly dierent from that of SARS-CoV-1 in terms of the number of subgraphs, particularly the bigger ones, that they have. Results of the F -test and KS-test for dierentiating the main proteases of SARS-CoV-1 and of SARS-CoV-2 based on Tr (E α,1 ) and Ῡ with two values of the fractional parameter α. Signicant values for p = 0.95 are boldfaced.

	test	α = 0.75 Tr (E α,1 )	Ῡ	α = 0.5 Tr (E α,1 )	Ῡ
	di. (%)	-48.66	-47.16	-2013.3	-1943.6
	F -test	0.6926	0.6843	0.9998	0.9998
	KS-test	0.9795	0.9795	0.9795	0.9795

  we resume of the studies reported for drug repurposing based on the main protease

	of SARS-CoV-2 as a target.
	template	drugs found	ref.
	N3-Mpro	cinanserin	[119]
	N3-Mpro	rendesivir, darunavir,	[128]
		saquinavir
	ketoamide-Mpro	sulfamethizole, sulfathiazole,	[234]
		kanamycin, tobramycin,
		phthalylsulfathiazole;
		droperidol, eszopiclone
		homotropine; alpelisib;
		tizanidine; mannitol
	N3-Mpro	Zanamivir, Indinavir,	[109]
		Saquinavir, and Remdesivir
	N3-Mpro	d-viniferin, myricitrin,	[121]
		taiwanhomoavone A,
		lactucopicrin 15-oxalate,
		nympholide A, afzelin,
		biorobin, hesperidin ,
		phyllaemblicin B

Table 10 :

 10

Some examples of drug repurposing experiments based on inhibiting the main protease of SARS-CoV-2.

Table 11 :

 11 Values of ∆G α ij , L and N BS based on two values of the fractional parameter

Table 15 :

 15 Resume of some of the compartmental models using network approaches described in the literature to analyze COVID-19 epidemiological parameters.

	Authors	Model	Findings	Ref.
	Prasse et al.	SIR + networks	evolution of epidemics with	[189]
			disaggregated β
	Small &	SEIR + networks	necessity of including contact	[219]
	Cavanagh		networks
	Linka et al.	SEIR +	good predictive capacity at earlier	[150]
		transport	stage of the epidemics
		networks
	Chung &	SEIR +	evolution of R (t) in Singapore	[49]
	Chew	multiplex
	Peirlinck et	SEIR + networks	reports latent and contact periods	[184]
	al.		in China and USA
	Block et al.	SEIR + networks	eects of reducing social contacts	[31]
			on epidemics progression
	Aleta and	SEIR +	epidemics progression at	[8]
	Moreno	metapopulation	region-level in Spain
		+ data analysis
	Chinazzi et	extended SEIR +	eects of traveling from China on	[48]
	al.	traveling network	the progress of the epidemics on
			the world

Table 17

 17 we report these values and in parenthesis the number of studies used for obtaining them. In general, it is seen that as average an infected individual infects 3 susceptible ones. This number is bigger than the ones reported for H1N1 inuenza in South Africa in 2009, Ebola in Guinea in 2014 and Zika in South America in 2015-16, but smaller than the ones reported for SARS epidemic in 002-03, smallpox in Indian subcontinent in 1968-73, poliomyelitis in Europe in 1955-60 and measles in Ghana in 1960-68 (see Table

Table 17 :

 17 Resume of the results of a meta-analysis for the basic reproduction number and case fatality rate (CFR) of COVID-19 in dierent regions of the world.ing studies carried out for COVID-19. Brauner et al.[START_REF] Brauner | Inferring the eectiveness of government interventions against COVID-19[END_REF] have estimated the eects of NPI on COVID-19 transmission in 41 countries during the rst wave of the pandemic.Their results are resumed in Fig.39where by using a Bayesian hierarchical model the authors linked NPI implementation dates to national case and death counts. As can be clearly seen closing all educational institutions, limiting gatherings to 10 people or less, and closing face-to-face businesses produced the biggest reduction of the eective reproductive number, here denoted by R t . However, the additional eect of stay-at-home orders was comparatively small.

	Gatherings limited to 1000 people or less
	Gatherings limited to 100 people or less
	Gatherings limited to 10 people or less
	Some businesses closed
	Most nonessen al businesses closed
	Schools and universi es closed
	Addi onal benefit of stay-at-home order
	on top of above NPIs
	0.0%	17.5%	35.0%	52.5%
	Posterior median reduc on in R t
	Fig. 39: Resume of the eectiveness of non-pharmaceutical interventions estimated
	across a suite of 206 analyses with dierent epidemiological parameters, data, and
	modeling assumptions as reported by Brauner et al. Reproduced from J. M. Brauner
	et al., Science 371, eabd9338 (2021) [37], published under CC BY 4.0.
	10 Non-pharmaceutical interventions with minimal interruptions

  (c). Therefore, a SIS model is an appropriate one for modeling this scenario.We should remark that the SIS model has been applied to similar situations for instance byOmi¢ and Van Mieghem [182] who considered a city as a node which can have infection introduced over the air transportation network from other cities. In this scenario they considered the SIS model as the adequate modeling tool. In another workMatamalas et al. [159] used again SIS for the worldwide air transportation network, with the goal of identifying the most important connections between airports for the spreading of epidemics and evaluate the epidemic incidence after its deactivation.Sanders et al. [205] analyzed the spread of an infection disease through 23

subpopulations via (documented) air trac data, and considering that the country is internationally quarantined. Similarly,

Qu and Wang [191]

,

Scaman et al. [207]

,

Onoue et al. [183]

, and

Ruan et al. [202]

,

Meloni et al. [162]

and

Ye et al. [256]

, among others, used SIS for modeling diseases propagating through airports in a network system. 10.1 Communicability function and the solution of the SIS model 10.1.1 Motivation

Here we mainly focus in nding the potential main routes of transmission of a viral disease through the nodes and edges of a network. That is, considering an epidemic disease which is transmitted through an airport network via a SIS mechanism we are interested in nding the shortest paths used by the epidemic to go from one airport to another. Due to the diusive nature of epidemic propagation it is not expected that such paths coincide in general with the shortest topological paths. Therefore we start here by considering a structural interpretation of the SIS model in terms of the communicability function, which can be directly linked to diusive paths in a network.

  That is, if Ŝv (t, ev) is relatively high it indicates that this node remains infected for some time, possibly due to continuous reinfections. This reinfection capacity can also be calculated for any node w as: Ŝw (t, ew) = e ζtA ew

w . Now, if we consider the term Ŝw (t, ev) = e ζtA ev w ,

Table 18 :

 18 Results of the edge-removal exercise for the airport connection networks of the UK.

	time	route	∆t ⋆ (%)	∆ C (%)	∆ l (%)
	200	Heathrow-Edinburgh	38.99	-9.93	0.049
	400	Heathrow-Manchester	88.89	-9.94	0.152
	600	Heathrow-Glasgow	194.44	-5.61	0.254
	1000	Heathrow-Dublin City	333.33	-3.60	0.407

Acknowledgements The author thanks nancial support from Ministerio de Ciencia, Innovacion y Universidades, Spain for the grant PID2019-107603GB-I00 Hubs-repelling/attracting Laplacian operators and related dynamics on graphs/networks.

Switch and restart mechanism

The subdiusive process emerging from the application of the FDE to study the propagation of perturbations through a host PPI network has a very small rate of convergence to the steady state. This creates the problem, as remarked by Batada et al. [START_REF] Batada | Stochastic model of proteinprotein interaction: Why signaling proteins need to be colocalized[END_REF], that such (sub)diusive processes alone are not sucient for carrying out cellular processes at a signicant rate in cells. However, this process produces dramatic increment of the probability that certain proteins are perturbed at very short times. This kind of shock wave eect of the transmission of perturbations occurs at much earlier times in the subdiusive regime than at the normal diusion one [START_REF] Estrada | Topological analysis of SARS CoV-2 main protease[END_REF]. Therefore, we have proposed a switch and restart process in which a subdiusive process starts at a given protein of the PPI, perturbs a few others, which then become the starting point of a new subdiusive process and so on as illustrated in Fig. 30. In Fig. 30 (a) a host protein (A) which has been previously perturbed by a viral protein navigates subdiusively inside the host cell until it nds host protein (B). If the response of B to its interaction with A is large enough in a short time, it can be considered as a perturbed protein. Then, the process is switched to start again at the protein (B). This is illustrated in (b) where (B) navigates subdiusively until it nds protein (C), which may also be perturbed. Fig. 30: Illustration of the switch-and-restart mechanism proposed here to explain the transmission of perturbations between proteins in a PPI (see text for explanation).