Julien Riposo

M Julien Mozziconacci

M Carlo Bianca

Anne Mme

Gégout-Petit

Rapporteur M Benjamin Audit

M Santo Motta

M Jean-Pierre

Indjehagopian Professeur

Examinateur M Jean-Bernard

Zuber

Computational and Mathematical Methods for Data Analysis in Biology and Finance

Introduction

Since the last century, the properties of the adjacency matrix M of a graph G (matrix whose element M ij is 1 if the nodes i and j are connected, and 0 otherwise) and the degree vector of the graph (vector whose element i is the number of connections of the node i with the other nodes of the graph) have been analyzed. Specifically the relationship between the principal eigenvector of M (the eigenvector related to the highest eigenvalue) and the degree vector of G has been investigated in, see among others, papers [START_REF] Mieghem | Double orthogonality and the nature of networks[END_REF][START_REF] Lovasz | Random walks of graphs: A survey[END_REF][START_REF] Serre | Matrices: Theory and Applications[END_REF] and constitutes one of the main interest of the theory of graph spectra [START_REF] Cvetkovic | Spectra of Graphs[END_REF][START_REF] Cvetkovic | An Introduction to the Theory of Graph Spectra[END_REF][START_REF] Germina | On products and line graphs of signed graphs, their eigenvalues and energy[END_REF]. In particular, it has been shown that it is possible to construct the original adjacency matrix from the principal eigenvector [START_REF] Porto | Reconstruction of protein structures from a vectorial representation[END_REF]. The relation between the degree vector of a graph G and the principal eigenvector of the related adjacency matrix M has been also proven to be of great interest in other research fields such as biology and finance. Indeed the spectrum of the adjacency matrix performs an important task considering the important role that has been shown recently in the dynamic of different phenomena occuring on biological and social systems, e.g. virus spread [START_REF] Pastor-Satorras | Epidemic dynamics and endemic states in complex networks[END_REF], vehicular traffic [START_REF] Darwish | raffic density estimation in vehicular ad hoc networks: A review[END_REF], social communication [START_REF] Coyle | Social networking: Communication revolution or evolution?[END_REF]. Specifically the vertices of the graph represent the individuals or the organizations and the edges of the graph mimic the interactions among them. In particular the principal eigenvector related to the largest eigenvalue (the spectral radius of the graph G) of the adjacency matrix of G characterizes dynamics process on networks. Indeed the ith-component of the principal eigenvector tends to be large if the ith-node has a large degree or if the components corresponding to its neighbors are large. Even if the topological meaning of the spectral radius is far from being well understood, it is nowadays recognized that the best strategy to minimize the spectral radius by links/nodes removal is based on the principal eigenvector [START_REF] Li | Bounds for the spectral radius of a graph when nodes are removed[END_REF]. Further applications refer to the componentwise-sorted principal eigenvector which contains information on the cluster topology and the 3D geometry of proteins [START_REF] Pietrucci | Graph theory meets ab initio molecular dynamics: Atomic structures and transformations at the nanoscale[END_REF]. Bearing all above in mind, looking at the relation between the principal eigenvector and the degree vector can be an approach for analyzing the data in biology and finance. This is the aim of the first part of the present thesis. A different approach for the treatment of data furnished by the financial or biological field is the time series analysis. In order to perform a good treatment of this data, different timediscrete models (difference equations) have been proposed. The most famous time-discrete model dates back to Robert Engle who has introduced in [START_REF] Engle | Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation[END_REF] the AutoRegressive Conditional Heteroskedasticity (ARCH) process for the study of various economic and financial data series of many countries. In this model the conditional variance (squared volatility) is not constant over time and shows autoregressive structure. Subsequently Bollerslev has proposed in [START_REF] Bollerslev | Generalized autoregressive conditional heteroscedasticity[END_REF] the Generalized ARCH (GARCH) process which provides a parsimonious parametrization for the conditional variance. Recently, an attempt to the definition of the corresponding timecontinuous limit for ARCH and GARCH has been proposed in [START_REF] Nelson | ARCH models as diffusion approximations[END_REF][START_REF] Kallsen | COGARCH as a continuous-time of GARCH(1, 1)[END_REF][START_REF] Chen | Equity volatility smile and skew under a CEV-based structural leverage model[END_REF] and this constitutes another interest of the thesis. The analysis of biological and financial systems can be also pursued by employing mathematical modeling by means of theoretical physics tools. The most famous approach is the ordinary differential equation-based model, in brief, ODE-based model. This approach assumes that the system (composed by human cells, molecules) is homogeneous thus ignoring the spatial structure of the system at the microscopic scale. In particular, ODE-based models have been employed to model cell growth, cell proliferation and the competition with the immune system, see the review [START_REF] Eftimie | Interactions between the immune system and cancer: a brief review of non-spatial mathematical models[END_REF] and the reference therein. Moreover, in order to take into account that some processes are CONTENTS not instantaneous but take a finite time to occur, a delay (or time-lag) can be considered, see [START_REF] Baker | Mathematical modelling of the interleukin-2 T-cell system: a comparative study of approaches based on ordinary and delay differential equations[END_REF] and therein references. Usually finding the analytical solution is difficult, then the analysis of the system refers to the existence of steady states and the related stability conditions [START_REF] Waniewski | Mathematical Modeling of the Immune Response[END_REF]. The modeling of a biological or financial system can be also described by means of suitable partial differential equations for the macroscopic variables (typically mass, momentum and energy). These equations are usually derived by following the guidelines oriented by continuum mechanics approach and specifically by employing mass, momentum, and energy conservation equations properly closed by phenomenological models corresponding to the material behavior of the system, see books [START_REF] Mollica | Modeling of Biological Materials[END_REF][START_REF] Chauviere | Cell Mechanics[END_REF], the review paper [START_REF] Fusi | Macroscopic models for pbroproliferative disorders -a review[END_REF] and the references cited therein. Recently, a new modeling approach has been proposed for the modeling of complex systems based on a generalized kinetic approach, called the thermostated kinetic theory [START_REF] Bianca | Thermostatted kinetic equations as models for complex systems in physics and life sciences[END_REF]. According to this approach the overall system is composed of particles, called active particles, whose microscopic state includes a scalar variable (activity) which models the function/strategy expressed by the particles. The particles expressing the same function/strategy are collected into subsystems, called functional subsystems. The description of each functional subsystem is based on a statistical approach by employing distribution functions on the particle microscopic state. Different types of interaction which have, as a result, activity transition, proliferation, destruction, and mutation, are responsible for the evolution of the system. This approach has been already employed for the modeling of biological systems, economic systems, social systems, politics and conflicts, see the reference section of the book [START_REF] Bianca | Towards a Mathematical Theory of Multiscale Complex Biological Systems[END_REF]. The modeling of a possible therapy for keloid formation by employing this approach is also treated in this thesis. This thesis is divided into two parts. The first part is concerned with a rigorous analysis on the relationship between the principal eigenvector and the degree vector of a matrix (positive, symmetric and irreducible) and its applications to the biological and financial fields. The second part is devoted to the derivation and analysis of mathematical models for wound healing diseases by means of methods from theoretical physics. Moreover the derivation of some properties for the continuous limit of the discrete stochastic processes (ARCH and GARCH) has also been collected in this part.

Specifically the contents of this thesis are organized as follows.

The chapter 1 collects the basic definitions and the main results on the matrix theory and the graph theory with the aim to be an introductory chapter for the reader. Specifically, the Perron-Frobenius theorem, the Sylvester's Inertia Law, the Gerschgoring theorem, the Semi-circle and the Markovitz-Pastur laws have been reviewed in this chapter. Moreover, the assortativity (a measure of the connections between high-degree nodes) is also introduced in this chapter. Finally the reader is led to the main problem of the thesis which consists of investigations of the behavior of the principal eigenvector. It is worth stressing that this chapter does not contain new results and the presentation does not aim at being exhaustive from the mathematical point of view. The interested reader in a more deeper understanding of the contents of this chapter will be referred to the pertinent literature.

The chapter 2 is mainly concerned with the derivation of a new upper bound for the error vector, which is the difference between the principal eigenvector of a matrix M (positive, symmetric and irreducible) and the degree vector of M (vector whose element i is the sum of the elements of the ith column of the matrix M ). Moreover a sufficient and necessary condition under which the principal eigenvector of M is proportional to the degree vector of M , is established. The previous analysis is generalized to the case of eventually non-negative matrices. Moreover, with the aim of depicting the analytical results, this chapter also contains some numerical simulations on different known graphs: the Erdos-Reyni graph (uniformly connected node) [START_REF] Erdos | On random graphs[END_REF], the Scalefree graph (power law degree distribution) [START_REF] Albert | Statistical mechanics of complex networks[END_REF], the Bipartite graph (the vertices set is divided into two subsets S 1 and S 2 with the property that there are no connections among the nodes of S 1 and among the nodes of S 2 ). A novel graph is defined in this chapter: the cherckerboard graph (a graph where the adjacency matrix is a checkerboard where the black boards are mainly composed of 1's and the white boards are mainly composed of 0's). Some applications to the field of biology, with genomic contact maps, and finance, with an autocorrelation matrix, are also shown in this chapter. The main results of this chapter have been collected in the papers [START_REF] Riposo | On the adjacency matrix of graphs: Principal eigenvector versus degree vector[END_REF] and [START_REF] Riposo | Principal eigenvector versus degree vector: Applications to networks[END_REF].

The chapter 3 deals with the statistical analysis of HiC matrices M , which are genomic contact maps where the element M ij is 1 if the ith and jth locus are in physical contact in the nucleus, and 0 otherwise. The aim of the statistical investigation is to find the nuclear compartments, an idea firstly introduced into the paper [START_REF] Lieberman-Aiden | Comprehensive mapping of long range interactions reveals folding principles of the humain genom[END_REF]. The chapter is organized as follows: first we introduce the HiC concept and focus on the outlines of the experimentation procedure, from the viewpoint the paper [START_REF] Lieberman-Aiden | Comprehensive mapping of long range interactions reveals folding principles of the humain genom[END_REF]. Secondly we give a sketch of ideas why the highest eigenvalue is high for HiC matrices, and then we focus on the statistical analysis of specific HiC matrices: one obtained from the human Embryonic Stem Cell (hESC) and the other one from the fibroblast. The statistical methods employed in this chapter are the kmeans on the 2D space of the second and third eigenvectors, with informational criteria. The repartition of biological features, as the Transcription Factors (TFs) or the histone marks, over the compartments is also analyzed. The main results of this chapter have been collected into the paper [START_REF] Lesne | 3D genome reconstruction from chromosomal contacts[END_REF].

The chapter 4 deals with further developments on a mathematical model recently proposed in [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] for the modeling of the acute inflammatory response to infection or trauma. In particular in order to take into account that some interactions have not an immediate effect, we introduce time delays. Specifically the chapter deals with the existence of steady states, determining the parameter regimes where the equilibrium points are stable, and the onset of Hopf bifurcation appears. Numerical simulations are performed with the main aim of supporting the analytical results and investigate further dynamics. It is worth pointing out that this chapter is an extension of paper [START_REF] Bianca | A delayed mathematical model for the acute inflammatory response to infection, to appear[END_REF].

The chapter 5 deals with the modeling of a wound healing disease under a therapeutic action by employing the methods of the thermostated kinetic theory for active particles [START_REF] Bianca | High-order moments conservation in thermostatted kinetic models[END_REF]. In particular, in order to test a therapeutic action for the keloid formation and the possible development of a cancer, an external force field coupled to a Gaussian thermostat is introduced into a mathematical model recently proposed in the paper [START_REF] Bianca | Mathematical modelling for keloid formation triggered by virus: Malignant effects and immune system competition[END_REF]. Specifically the model depicts the competition of the immune system cells with a virus, the mutated fibroblast cells, and the cancer cells. Employing a computational analysis, the effect of three different external force fields mimic therapeutic actions is analyzed: A vaccine for the virus, the PUVA therapy for the keloid and a vaccine for the cancer. The results are in agreement with the evidence that the sole action of the immune system is not sufficient to obtain a total depletion of keloid. Further refinements and developments are also discussed into the chapter. The contents of this chapter can be found in the paper [START_REF] Bianca | Mimic therapeutic actions against keloid by thermostatted kinetic theory methods, to appear[END_REF].

Part I

Methods for Data Analysis in

Biology and Finance

Chapter 1

The Mathematical Tools

This chapter deals with the main definitions and results on matrix theory and graph theory, and specifically on the random matrix, sampled matrix and the assortativity. The presentation does not aim at being exhaustive from the mathematical point of view. Indeed the definitions are limited to the contents of the present thesis. The interested reader in a deeper understanding on contents of this chapter will be referred to the pertinent literature.

Backgrounds on matrix algebra

Let K be a scalar field, n, m ∈ N * = {a ∈ N, a > 0}, M n,m (K) the vector space of the matrices M = (M ij ) 1 i n,1 j m of dimension n × m whose elements are such that M ij ∈ K, and M n (K) the vectorial space of the square matrices of dimension n 2 .

Definition 1.1.1 A matrix M ∈ M n (R) is said • Nonnegative, or positive, if M ij 0, ∀i, j ∈ {1, 2, . . . , n};

• Reducible if there exists a permutation matrix P such that

P M t P = A 11 A 12 0 A 22 ,
where A 11 ∈ M r (R), M 22 ∈ M n-r (R), and A 12 ∈ M r,n-r (R), with 0 < r < n.

A matrix M is said to be irreducible if it is not reducible.

In what follows, we denote by SIR n (K) the set of symmetric and irreducible matrices, whose elements are such that M ij ∈ K, ∀i, j ∈ {1, 2, . . . , n}, by R + = {a ∈ R, a 0}, R * + = {a ∈ R, a > 0}, t M the transpose of M , and t e = (1, 1, . . . , 1) ∈ R n . Note that M e ∈ R n is the vector whose generic element (M e) (i) is the sum of the elements belonging to the ith-row of M , namely

(M e) (i) = n j=1 M ij .
Let M ∈ SIR n (R + ). Since M is a real symmetric matrix, by the spectral theorem [START_REF] Cauchy | Sur l' équation à l'aide de laquelle on détermine les inégalités séculaires du mouvement des planètes[END_REF], there exists an orthonormal basis B M = {V 1 [M ], V 2 [M ], . . . , V n [M ]} of eigenvectors of M , related to the real eigenvalues {Λ 1 [M ], Λ 2 [M ], . . . , Λ n [M ]}. In particular, we assume that

|Λ 1 [M ]| |Λ 2 [M ]| • • • |Λ n [M ]|.
(1. [START_REF] Mieghem | Double orthogonality and the nature of networks[END_REF] In what follows, we denote by λ i [M ], for i ∈ {1, 2, . . . , n}, the eigenvalues of the matrix M arranged as follows:

λ 1 [M ] λ 2 [M ] • • • λ n [M ]. (1.
2)

The basis

B M = {V 1 [M ], V 2 [M ], . . . , V n [M ]} is thus such that < V i [M ], V j [M ] >= 1 if i = j 0 if i = j (1.3)
where < • , • > denotes the canonical inner product in R n . In particular

V i [M ] = (V (1)
i [M ], V (2)
i [M ], . . . , V

(n) i [M ]) is such that < V i [M ], V i [M ] >= n j=1 (V (j) i [M ]) 2 = 1.
(1.4)

Note that

< V i [M ], e >= n j=1 V (j) i [M ].
The following result fulfills an important role in the present thesis.

Theorem 1.1.1 (Perron-Frobenius [START_REF] Perron | Zur Theorie der Matrices[END_REF], [START_REF] Frobenius | Uber matrizen aus positiven elementen[END_REF]) Let M ∈ SIR n (R + ) and ρ(M ) the spectral radius of M given by ρ(M ) = max(|Λ| : Λ ∈ {Λ 1 [M ], ..., Λ 2 [M ]}).

(1.5)

Then

• ρ(M ) > 0 is the largest eigenvalue of M , i.e. ρ(M ) = λ 1 [M ] > 0;

• The algebraic multiplicity of ρ(M ) is 1;

• The components of the eigenvector associated with ρ(M ) are strictly positive;

• The only eigenvector of M that has strictly positive components is the eigenvector associated with ρ(M ).

It is worth noting that λ

1 [M ] = |Λ 1 [M ]|.
Bearing all above in mind and since M ∈ SIR n (R + ), by the Perron-Frobenius theorem 1.1.1 we have

Λ 1 [M ] > 0, V (j)
1 [M ] > 0, ∀j ∈ {1, 2, . . . , n}.

(1.6)

In particular the algebraic dimension of Λ 1 [M ] is 1. Accordingly Eq (1.1) now reads:

Λ 1 [M ] |Λ 2 [M ]| • • • |Λ n [M ]|, (1.7)
and Eq (1.

2) now reads:

λ 1 [M ] > λ 2 [M ] • • • λ n [M ]. (1.8) Remark 1.1.1 The algebraic multiplicity of Λ 1 [M ] is 1. Therefore Λ 1 [M ] = |Λ 2 [M ]| implies Λ 2 [M ] = -Λ 1 [M ]
< 0 and then by Eq (1.8) one has

λ n [M ] = Λ 2 [M ].
It is well-known that (see for instance [START_REF] Meyer | Matrix analysis and applied linear algebra[END_REF]) min i∈{1,2,...,n} (M e) (i) ρ(M ) max i∈{1,2,...,n} (M e) (i) .

(1.9) (M e) (i) (M is called regular).

Example 1.1. [START_REF] Mieghem | Double orthogonality and the nature of networks[END_REF] The following matrix M ∈ M 2 (R)

M = 0 1 1 1
is a non-negative and symmetric matrix. It is also irreducible, indeed the only permutation matrices are P 1 = 1 0 0 1 = Id and P 2 = 0 1 1 0

and

P 1 M t P 1 = M, P 2 M t P 2 = 1 1 1 0 . Moreover min i∈{1,2}
(M e) (i) = 1 and max i∈{1,2}

(M e) (i) = 2. Then α = 1 and β = 2, and 1 (M e) (i) 2 ∀i ∈ {1, 2}.

Definition 1.1.3 Let M ∈ SIR n (R + ). We define

m = 1 αβ |Λ 2 [M ]| Λ 1 [M ] , (1.12)
where α and β are given by Eq (1.10) and Eq (1.11), respectively.

Theorem 1.1.2 (Sylvester's Inertia Law [START_REF] Sylvester | A demonstration of the theorem that every homogeneous quadratic polynomial is reducible by real orthogonal substitutions to the form of a sum of positive and negative squares[END_REF]) Let A, P ∈ M n (R) with P invertible. Then the matrix A has the same numbers of positive, negative and zero eigenvalues as t P AP .

Under the assumption of the theorem 1.1.2, if P is also a symmetric matrix then t P AP = P AP = P -1 (P 2 A)P . In particular the matrices P 2 A and P -1 (P 2 A)P have the same eigenvalues. Moreover the matrices P 2 A and P AP have the same eigenvalues. Then the matrices A and P 2 A have the same number of positive, negative, and zero eigenvalues.

Let M ∈ SIR n (R + ) and D ∈ M n (R + ) the following diagonal matrix whose elements are

D ij =      n k=1 M ik if i = j 0 if i = j (1.13)
Obviously D is a symmetric invertible matrix and the eigenvalues of D read:

λ i [D] = n k=1
M ik > 0, ∀i ∈ {1, 2, . . . , n}.

Proposition 1.1.1 Let M ∈ SIR n (R + ). Then 1) D is a positive-definite matrix; 2) The matrix M = D -1 M, (1.14)
where

D -1 ij =            1 n k=1 M ik if i = j 0 if i = j (1.15)
is well-defined. In particular M is such that

M ij = M ij (M e) (i) = M ij n k=1 M ik , (1.16
)

and then 0 M ij 1, n k=1 M ik = 1 (M is called stochastic matrix).
3) The ith-eigenvalue of D -1 reads:

λ i [D -1 ] = 1 n j=1 M ij , ∀i ∈ {1, 2, . . . , n} .
4) The characteristic polynomial of the matrix D 1/2 M D -1/2 is the same of the matrix M , where

D -1 2 ij =            1 n k=1 M ik if i = j 0 if i = j (1.17)
5) The matrices M and M have the same number of positive, negative and zero eigenvalues.

Proposition 1.1.2 Let M ∈ SIR n (R + ) and M = D -1 M the matrix defined in (1.14). Then M is irreducible.

Proof Suppose that M is a reducible matrix, then there exists a permutation matrix P such that

P M t P = A 11 A 12 0 A 22 ,
where A 11 ∈ M r,r (R), A 22 ∈ M n-r,n-r (R), and A 12 ∈ M r,n-r (R), with 0 < r < n. The null matrix has dimension (n -r) × r.

Since M = DM we have M ij = D ii M ij , ∀i, j ∈ {1, 2, . . . , n}.
It turns out that we also have

P M t P = A 11 A 12 0 A 22 ,
where A 11 ∈ M r,r (R), A 22 ∈ M n-r,n-r (R), and A 12 ∈ M r,n-r (R), with 0 < r < n. Then M is reducible.

Theorem 1.1.3 (Courant-Fischer or Min-Max, [START_REF] Reed | Methods of modern mathematical physics iv: Analysis of operators[END_REF]) Let M ∈ M n (R) be a symmetric matrix and S ⊂ R n a proper subspace of R n . Then

λ i [M ] = max dim(S)=i min x =0 x∈S < M x, x > < x, x > , ∀i ∈ {1, 2, . . . , n}.
In particular

λ 1 [M ] = max x =0 < M x, x > < x, x > and λ n [M ] = min x =0 < M x, x > < x, x > .
Theorem 1.1.4 (Gerschgorin [START_REF] Gerschgorin | Uber die Abgrenzung der Eigenwerte einer Matrix[END_REF]) Let M ∈ M n (C). Then for any eigenvalue λ ∈ C of M , there exists i(λ) ∈ {1, 2, . . . , n} such that

|λ -M i(λ)i(λ) | n j=1,j =i(λ) |M i(λ)j |. (1.18)
Remark 1.1.4 By the triangular inequality, the Eq (1.18) reads:

|λ| -|M i(λ)i(λ) | |λ -M i(λ)i(λ) | n j=1,j =i(λ) |M i(λ)j |,
and then

|λ| n j=1 |M i(λ)j |.

Background on graph theory

This section is devoted to the main definitions and properties about graph theory. The following result holds.

I = {i 1 , i 2 , . . . , i s }, J = {j 1 , j 2 , . . . , j t }, s + t = n, ∀s, t 1, such that, ∀(i, j) ∈ I × J, M ij = 0. A matrix M is irreducible if it is not reducible.
Proof Assume that there exists a partition of the set {1, 2, . . . , n} into two subsets:

I = {i 1 , i 2 , . . . , i s }, J = {j 1 , j 2 , . . . , j t }, s + t = n, ∀s, t 1,
such that, ∀(i, j) ∈ I × J, M ij = 0. Let t P be the permutation matrix which orders the columns of M in the order (j 1 , . . . , j t , i 1 , . . . , i s ):

M =    M 11 • • • M 1n . . . . . . . . . M n1 • • • M nn    , then
M t P =    M 1j 1 • • • M 1jt M 1i 1 • • • M 1is . . . . . . . . . . . . . . . . . . M 1j 1 • • • M 1jt M 1i 1 • • • M 1is    .
In order to apply the same permutation on the rows of the matrix M , we multiply on the left of M t P by P :

P M t P =           M j 1 j 1 • • • M j 1 jt M j 1 i 1 • • • M j 1 is . . . . . . . . . . . . . . . . . . M jtj 1 • • • M jtjt M 1i 1 • • • M jtis M i 1 j 1 • • • M i 1 jt M i 1 i 1 • • • M i 1 is . . . . . . . . . . . . . . . . . . M isj 1 • • • M isjt M isi 1 • • • M isis           P M t P =           M j 1 j 1 • • • M j 1 jt M j 1 i 1 • • • M j 1 is . . . . . . . . . . . . . . . . . . M jtj 1 • • • M jtjt M 1i 1 • • • M jtis 0 • • • 0 M i 1 i 1 • • • M i 1 is . . . . . . . . . . . . . . . . . . 0 • • • 0 M isi 1 • • • M isis           .
This last matrix has the seeked shape. Conversely, there exists a permutation matrix P such that

P M t P = A 11 A 12 0 A 22 ,
where

A 11 ∈ M s (R), A 22 ∈ M t (R),
A graph G = (V, E) is said connected if ∀u, v ∈ V,
there exists a path from u to v in G.

Definition 1.2.5 (Adjacency matrix M of a graph) Let G = (V, E) be a graph of order n.

The adjacency matrix of G is the n × n symmetric matrix M = (M ij ) 1 i,j n with entries

M ij = 1 v i v j ∈ E 0 v i v j / ∈ E Proposition 1.2.2 A graph G = (V, E) is connected if and only if its adjacency matrix M ∈ SIR n (R + ).
Proof First we show that if M is reducible then G is not connected. Suppose M reducible.

Then there exists a permutation matrix P ∈ M n (R) such that:

P M t P = A 11 A 12 0 A 22 ,
where

A 11 ∈ M s (R), A 22 ∈ M t (R)
, and A 12 ∈ M s,t (R), with s, t ∈ {1, 2, . . . , n} (especially the elements of the matrix A 12 are zeros). For i ∈ {t + 1, t + 2, . . . , n} and j ∈ {1, 2, . . . , t}, there is no edge between i and j in the graph whose adjacency matrix is P M t P . Then this graph is not connected, and G is not connected as well.

Conversely, Suppose that G is not connected. {1, 2, . . . , n} ∈ V is the set of vertices of G. There exists at least one couple (i 1 , i n ) / ∈ E such that there exists no path from i 1 to i n . Let

V 1 = {i 1 , i 2 , . . . , i t }
be the set of vertices of G such that there exists no path from i n to any of the elements of V 1 .

The other vertices are grouped in the set

V 2 = {i t+1 , i t+2 , . . . , i n-1 }.
By construction, one has:

∀(i, j) ∈ V 1 × V 2 , M ij = 0.
Let P be the permutation matrix corresponding to σ : (1, 2, . . . , n) → (i 1 , i 2 , . . . , i t , i t+1 , . . . , i n ).

Then P M t P is given by:

P M t P = A 11 A 12 0 A 22 ,
where

A 11 ∈ M s (R), A 22 ∈ M t (
M =     0 1 1 1 1 0 0 1 1 0 0 0 1 1 0 0     .
In particular M is irreducible.

Remark 1.2.2

The adjacency matrix M of a connected graph is such that M ∈ SIR n (R + ). In particular we have (see definition 1.1.2)

d max = αn, d min = n β .
(1.20)

Random matrices and sampled matrices

The reader interested in the proof of the results of this section is referred to [START_REF] Pastur | Eigenvalue Distribution of Large Random Matrices[END_REF]. Let (Ω, P) be a probability space.

Definition 1.3.1 A random variable X : Ω → K is a measurable function from the set of possible outcomes Ω to the scalar field K. The set of random variables from Ω to K is denoted

K Ω . Definition 1.3.2 A random matrix X = (X ij ) 1 i n,1 j m ,
where n, m ∈ N * , is a matrix whose elements are the random variables

X ij ∈ K Ω . The set of random matrices of dimension n × m is denoted by M n,m (K Ω ). Definition 1.3.3 The subset GOE n ⊂ M n (K Ω ) of symmetric matrices which elements are Gaussian random variables is called Gaussian Orthogonal Ensemble. Moreover, if X ∈ GOE n , then the distributions of the elements of X is invariant under orthogonal conjugation, that is, if X ∈ GOE n , O ∈ M n (K Ω ) such that t OO = Id n
, where Id n is the identity matrix of size n, then the distribution of the elements of X and of t OX O are the same.

Remark 1.3.1 Let B R be the sigma-algebra of the Borelians of R. For K = R, the eigenvalues of X have their values in a Borelian B ∈ B R (the set of Borelians of R is engendered by the intervals of the form ]a, +∞[, where a ∈ R).

Definition 1.3.4 Let X ∈ GOE n . The normalized counting measure N n associated to X is the number N n (B) = 1 n card{l ∈ {1, 2, . . . , n}, λ l [X ] ∈ B} ∀B ∈ B R . (1.21)
In what follows we denote by E[X ij ] the mathematical expectation of the random variable X ij :

E[X ij ] = Ω X ij (ω)dP(ω). (1.22) Definition 1.3.5 Let X = (X ij ) 1 i n,1 j m ∈ M n,m (K Ω ) such that E [X ij ] = 0, ∀i ∈ {1, 2, . . . , n}, ∀j ∈ {1, 2, . . . , m}. (1.23)
The covariance matrix C of X is the following n × n matrix:

C = 1 n X t X . (1.24) Theorem 1.3.1 (Semi-circle law) Let X = (X ij ) 1 i,j n ∈ GOE n be such that E [X ij ] = 0, E X 2 ij = (1 + δ ij )w 2 , ∀i, j ∈ {1, 2, . . . , n}, (1.25)
where w > 0 and δ ij is the delta of Kronecker. If N n is the normalized counting measure associated to X / √ n, then there exists a probability measure P such that (the convergence refers to n → +∞)

N n (B) P-a.s. -→ B ∆(λ)dλ ∀B ∈ B R , (1.26)
where

∆(λ) =        √ 4w 2 -λ 2 2πw 2 if |λ| 2w 0 otherwise (1.27)
The distribution given by Eq (1.27) has the shape of a semi-circle, hence the name.

Theorem 1.3.2 (Marchenco-Pastur law) Let X = (X ij ) 1 i n,1 j m ∈ M n,m (K Ω ) be such that E [X ij ] = 0, E X 2 ij = w 2 , ∀i ∈ {1, 2, . . . , n}, ∀j ∈ {1, 2, . . . , m}, (1.28)
where w > 0. Let Q = n/m, N n be the counting measure associated to C and

λ ± = w 2 1 + 1 Q ± 2 1 Q .
If N n is the normalized counting measure associated to C, then there exists a probability measure P such that N n (B)

P-a.s. -→ B ∆(λ)dλ ∀B ∈ B R , (1.29)
where ∆(λ) is the following Wishart distribution:

∆(λ) =        max 1 - 1 Q , 0 δ(λ) + Q 2πλw 2 (λ + -λ)(λ -λ -) if λ -λ λ + 0 otherwise (1.30)
where δ(λ) is the Dirac mass at 0.

The convergence is related to n → +∞ and m → +∞, but with

Q ∈ Q + . Remark 1.3.2 If X ∈ M n (K Ω ) then Q = 1 and ∆(λ) =        1 2πw 2 4w 2 -λ λ if 0 λ 4w 0 otherwise (1.31)
The functions defined in (1.27) and (1.30), called density of eigenvalues, represent the probability density that an eigenvalue of X is λ.

Bearing all of this in mind, we can define the sampled matrix associated to a random matrix.

Definition 1.3.6 Let X = (X ij ) 1 i n,1 j m ∈ M n,m (K Ω ) and ω ∈ Ω. The sampled matrix X s = (X s ij ) 1 i n,1 j m ∈ M n,m (K) related to ω is the matrix such that X s ij = X ij (ω). We write X s = X (ω).
Remark 1.3.3 From the n × n sampled matrix X s , we can define the empirical density of eigenvalues ∆ s (λ). This is actually the histogram of eigenvalues associated to the matrix X s . Remark 1.3.4 Any matrix M ∈ M n,m (K) can be associated to a random matrix X ∈ M n,m (K Ω ) with an existing ω ∈ Ω such that M = X (ω).

In the following we will denote by M , instead of X s , a sampled matrix when there will be not ambiguity.

The assortativity of a matrix

Definition 1.4.1 Let X ∈ SIR n (R Ω + ).
The quantity a D ((X e) (i) , (X e) (j) ) (D refers to 'Degree') related to the degrees (X e) (i) and (X e) (j) , ∀i, j ∈ {1, 2, . . . , n} is the correlation coefficient Cor[(X e) (i) , (X e) (j) ] between (X e) (i) and (X e) (j) :

a D ((X e) (i) , (X e) (j) ) := Cor[(X e) (i) , (X e) (j) ] = E[(X e) (i) (X e) (j) ] -E[(X e) (i) ]E[(X e) (j) ]
Var[(X e) (i) ] Var[(X e) (j) ] ,

(1.32) where Var[(X e) (i) 2 is the variance of the random variable (X e) (i) . a D is called local assortativity if the probability measure P takes into account the probability law of the graph connectedness. An estimator of the correlation coefficient Cor[(X e) (i) , (X e) (j) ] is the empirical correlation coefficient Cor[(M e) (i) , (M e) (j) ], whose empirical expectation must take into account the connectedness of the network.

] = E[((X e) (i) ) 2 ]-(E[(X e) (i) ])
Definition 1.4.2 ([45]) Let X ∈ SIR n (R Ω + ), ω ∈ Ω and M = X (ω) ∈ SIR n (R + ). The global assortativity A D (M ) of M is defined by A D (M ) = Cor[(M e) (i) , (M e) (j) ] = E[(M e) (i) (M e) (j) ] -(E[(M e) (i) ]) 2 E[((M e) (i) ) 2 ] -(E[(M e) (i) ]) 2 (1.33)
In particular we have

E[(M e) (i) ] = n i=1 n j=1 (M e) (i) M ij n i=1 n j=1 M ij = < M 2 e, e > < M e, e > , E[(M e) (i) (M e) (j) ] = n i=1 n j=1 (M e) (i) (M e) (j) M ij n i=1 n j=1 M ij = < M 3 e, e > < M e, e > , E[((M e) (i) ) 2 ] = n i=1 n j=1 ((M e) (i) ) 2 M ij n i=1 n j=1 M ij = n i=1 ((M e) (i) ) 3
< M e, e > .

The terms in the mean are weighted by the elements M ij , in order to take into account the connectivity of the graph. Replacing the above relations into Eq (1.33) yields

A D (M ) = A(M ) (< M e, e >) n i=1 ((M e) (i) ) 3 -(< M 2 e, e >) 2 , ( 1.34)
where

A(M ) = (< M e, e >)(< M 3 e, e >) -(< M 2 e, e >) 2 . (1.35)
Since the denominator of A D (M ) is associated to a variance one has < M e, e > n i=1 ((M e) (i) ) 3 -(< M 2 e, e >) 2 0. In particular the equality holds if M is a regular matrix. When the denominator is not zero, the sign of A D (M ) only depends on the sign of A, which will be the object of further investigations.

From the graph theory viewpoint when the assortativity is high then high-degree nodes are mainly connected with high-degree nodes [START_REF] Mieghem | Influence of assortativity and degree-preserving rewiring on the spectra of networks[END_REF] [START_REF] Li | Degree and principal eigenvectors in complex networks[END_REF].

The asymptotic behavior of

Λ 1 [M ] and V 1 [M ]
In this section we show some properties on sampled matrices introduced in the previous section: the semi-circle law Eq (1.27) and the Wishart law (Eq 1.30). The main aim is to introduce the reader to the main problem that will be debated in the first part of this thesis. In particular we will consider matrices whose size n ∈ {10, 50, 100, 150, 500, 1000} and whose elements are realizations of Gaussian random variables with expectation 0 and variance 1. Specifically we plot the empirical histogram of the related eigenvalues. According to Eq (1.27), Figure 1.2 shows that when n increases, the density curve (in red) approches a semi-circle. A possible explanation of these behaviors remains on the fact that when n increases the sampling is larger and the density approaches to the true value (this is in agreement with the central limit theorem for the density of the eigenvalues). Figure 1.3 shows the density of the eigenvalues related to the covariance matrix, obtained from each of the previous defined matrices. As Figure 1.3 shows, when n increases the density of the eigenvalues approaches to the Wishart distribution. Now we consider matrices whose size n ∈ {150, 500, 1000} and whose elements are realizations of Gaussian random variables with expectation 10 and variance 1. As Figure 1.4 shows, the magnitude of the eigenvalue Λ 1 [M ] increases and the difference between Λ 1 [M ] and Λ 2 [M ] increases when n increases. This is in agreement with the following lower and upper bounds [START_REF] Meyer | Matrix analysis and applied linear algebra[END_REF]:

min i∈{1,2,...,n}   n j=1 M ij   Λ 1 [M ] max i∈{1,2,...,n}   n j=1 M ij   . (1.36)
The aim of this part of the thesis is to obtain a new and more suitable relation between Λ 1 [M ] and the degree vector M e. Roughly speaking we decompose e = t (1, 1, . . . , 1) on the basis The red curves are related to Eq (1.31). According to Eq (1.31), the Wishart distribution is an asymptotical behavior.

B M = (V 1 [M ], V 2 [M ], . . . , V n [M ]) as follows: e = n i=1 c i V i [M ], where c i =< V i [M ], e >= n j=1 V (j) i [M ],
then we should obtain the following relation:

M e c 1 Λ 1 [M ]V 1 [M ].
(1.37)

From Eq (1.37) we have

< M e, M e > (c 1 Λ 1 [M ]) 2 < V 1 [M ], V 1 [M ] >= (c 1 Λ 1 [M ]) 2 ,
and then Bearing all above in mind, we have .38) According to the equation (1.38), the highest eigenvector V 1 [M ] is proportional to the degree vector M e, where the proportionality factor is 1/ √ < M e, M e >. We have seen in Figure 1.4 that when the mean of the elements is not zero but 10, Eq (1.37) holds true. In chapter 2, we will attempt to find a rigorous condition where Eq (1.38) is satisfied.

c 1 Λ 1 [M ] < M e, M e >.
V 1 [M ] 1 √ < M e, M e > M e. ( 1

Chapter 2 Adjacency Matrix: Principal and Degree Vectors

This chapter is concerned with the analysis of the difference between the principal eigenvector V 1 [M ] and the degree vector M e, where M ∈ SIR n (R + ) or M ∈ M n (R + ). Specifically we derive upper bounds for some matrix norms. Moreover we state a sufficient and necessary condition under which Eq (1.38) holds true. Finally some numerical simulations are performed on different known graphs, and applications to the field of biology and finance are shown. It is worth pointing out that this chapter is an extension of papers [START_REF] Riposo | On the adjacency matrix of graphs: Principal eigenvector versus degree vector[END_REF] and [START_REF] Riposo | Principal eigenvector versus degree vector: Applications to networks[END_REF].

Estimating the difference between the principal eigenvector and the degree vector

Definition 2.1.1 Let v = (v (1) , v (2) , . . . , v (n) ) be a vector of R n . The p th -order moment of v is defined as follows:

E p (v) = 1 n n j=1 (v (j) ) p . (2.1)
In particular, we are interested in the 1 st -order and 2 nd -order moments.

Bearing chapter 1 in mind and in order analytically estimate the commited error made by identfying the degree vector with the principal eigenvector, we consider the components of the vector e with respect to the orthonormal basis

B M = {V 1 [M ], V 2 [M ], . . . , V n [M ]}. Accordingly we have e = n i=1 c i V i [M ] = (c 1 , c 2 , . . . , c n ) B M . (2.2)
In particular we have

n =< e, e >=< n i=1 c i V i [M ], n i=1 c i V i [M ] >= n i=1 c 2 i . (2.3)
Moreover

M e = n i=1 c i M V i [M ] = n i=1 c i Λ i [M ]V i [M ]. (2.4)
Note that Eq (1.4) implies 0 < V

(j) < V 1 [M ], e >= n j=1 V (j) 1 [M ] > n j=1 (V (j) 1 [M ]) 2 = 1. (2.5) Since n j=1 V (j) i [M ] =< V i [M ], e >=< V i [M ], n j=1 c j V j [M ] >= c i , ∀i ∈ {1, 2, . . . , n}. (2.6)
In particular we have c 1 > 1. Then from Eq (2.4) we have

c 1 Λ 1 [M ]V 1 [M ] = M e - n i=2 c i Λ i [M ]V i [M ],
and since Λ 1 [M ] > 0 and c 1 > 1 > 0, we have

1 c 1 Λ 1 [M ] n i=2 c i Λ i [M ]V i [M ] = M e c 1 Λ 1 [M ] -V 1 [M ]. (2.7)
Definition 2.1.2 We define the error vector E as follows:

E = 1 c 1 Λ 1 [M ] n i=2 c i Λ i [M ]V i [M ].
(2.8)

The following result holds true.

Theorem 2.1.1 Let M ∈ SIR n (R + ), α and β defined by Eq (1.10) and Eq (1.11), and m defined by Eq (1.12). If m(αβ) 2 < 1 then

|E 1 (E)| m √ n (αβ -1) αβ (1 -m(αβ) 2 )(1 -mαβ) E 2 (E) m √ n αβ αβ -1 1 -m(αβ) 2 (2.9) Proof Since |Λ i [M ]| Λ 1 [M ] |Λ 2 [M ]| Λ 1 [M ] , ∀i ∈ {2, 3, . . . , n}, (2.10
)

and c i =< V i [M ], e >= n j=1 V (j) i [M ], ∀i ∈ {1, 2, . . . , n}, one has |E 1 (E)| = 1 n 1 c 1 Λ 1 [M ] n j=1 n i=2 c i Λ i [M ]V (j) i [M ] = 1 c 1 n n i=2 c i Λ i [M ] Λ 1 [M ]   n j=1 V (j) i [M ]   c i = 1 c 1 n n i=2 c 2 i Λ i [M ] Λ 1 [M ] 1 c 1 n n i=2 c 2 i Λ i [M ] Λ 1 [M ] (2.10) 1 c 1 n n i=2 c 2 i Λ 2 [M ] Λ 1 [M ] (2.3) = 1 c 1 n (n -c 2 1 ) Λ 2 [M ] Λ 1 [M ] .
Then

|E 1 (E)| 1 - c 2 1 n c 1 Λ 2 [M ] Λ 1 [M ] . (2.11)
By the same arguments we obtain

E 2 (E) 1 - c 2 1 n c 1 Λ 2 [M ] Λ 1 [M ] . (2.12) Since m(αβ) 2 < 1, we have mαβ < 1 αβ < 1.
In particular

Λ 2 [M ] Λ 1 [M ] = |Λ 2 [M ]| Λ 1 [M ] = mαβ < 1 αβ .

CHAPTER 2. ADJACENCY MATRIX: PRINCIPAL AND DEGREE VECTORS

Then n min i∈{1,2,...,n}

(M e) (i) n j=1 (M e) (j) = n j=1 n i=1 c i Λ i [M ]V (j) i [M ] = n i=1 c i Λ i [M ]   n j=1 V (j) i [M ]   = n i=1 c 2 i Λ i [M ] n i=1 |c 2 i Λ i [M ]| = n i=1 c 2 i |Λ i [M ]| =c 2 1 Λ 1 [M ] + c 2 2 |Λ 2 [M ]| + n i=3 c 2 i |Λ i [M ]| c 2 1 Λ 1 [M ] + c 2 2 |Λ 2 [M ]| + |Λ 2 [M ]| n i=3 c 2 i =Λ 1 [M ] c 2 1 + |Λ 2 [M ]| Λ 1 [M ] n i=2 c 2 i Λ 1 [M ] c 2 1 + mαβ n i=2 c 2 i =Λ 1 [M ] c 2 1 + mαβ(n -c 2 1
) . Therefore min i∈{1,2,...,n}

(M e) (i) Λ 1 [M ] c 2 1 n + mαβ 1 - c 2 1 n .
From Eq (1.9) we have

1 αβ = min i∈{1,2,...,n} (M e) (i) max i∈{1,2,...,n} (M e) (i) c 2 1 n + 1 - c 2 1 n mαβ,
and then

1 c 1 αβ n 1 -mαβ 1 -m(αβ) 2 . Moreover 1 - c 2 1 n 1 - 1 αβ 1 -m(αβ) 2 1 -mαβ ⇐⇒ 1 - c 2 1 n αβ -1 αβ(1 -mαβ) . ( 2

.13)

Bearing all above in mind, the inequalities (2.11) and (2.12) now reads:

|E 1 (E)| αβ n 1 -mαβ 1 -m(αβ) 2 αβ -1 αβ(1 -mαβ) mαβ = m √ n (αβ -1) αβ (1 -m(αβ) 2 )(1 -mαβ) , E 2 (E) αβ n 1 -mαβ 1 -m(αβ) 2 αβ -1 αβ(1 -mαβ) mαβ = m √ n αβ αβ -1 1 -m(αβ) 2 .
The proof is now concluded.

Remark 2.1.1 Under the assumption m(αβ) 2 < 1, it also holds

|Λ 2 [M ]| Λ 1 [M ] < min i∈{1,2,...,n} (M e) (i) max i∈{1,2,...,n}
(M e) (i) .

We now analyze the error (2.8) componentwise. From Eq (2.7), one has

V (j) 1 [M ] = (M e) (j) c 1 Λ 1 [M ] 1 - 1 (M e) (j) n i=2 c i Λ i [M ]V (j) i [M ] , ∀j ∈ {1, 2, . . . , n}, then 1 (M e) (j) n i=2 c i Λ i [M ]V (j) i [M ] = 1 -c 1 Λ 1 [M ] V (j) 1 [M ]
(M e) (j) , ∀j ∈ {1, 2, . . . , n}.

Definition 2.1.3 We define the error vector Ẽ as follows:

Ẽ(j) = 1 (M e) (j) n i=2 c i Λ i [M ]V (j) i [M ], ∀j ∈ {1, 2, . . . , n}. (2.14)
It is worth stressing that

| Ẽ(j) | β n n i=2 c i Λ i [M ]V (j) i [M ]
, ∀j ∈ {1, 2, . . . , n}.

Theorem 2.1.2 Let M ∈ SIR n (R + ), α and β defined by Eq (1.10) and Eq (1.11), and m defined by Eq (1.12). If m(αβ) 2 < 1, then

|E 1 ( Ẽ)| αβ -1 1 -mαβ mαβ E 2 ( Ẽ) αβ(αβ -1) 1 -mαβ mαβ (2.15)
Proof By using the same arguments of the proof of theorem 2.1.1, we have

|E 1 ( Ẽ)| 1 n β n n j=1 n i=2 c i Λ i [M ]V (j) i [M ] = β n n i=2 c 2 i Λ i [M ] n β|Λ 2 [M ]| n 1 - c 2 1 n .
One also has

E 2 ( Ẽ) β|Λ 2 [M ]| n 1 - c 2 1 n . Since Λ 1 [M ] max i∈{1,2,...,n} (M e) (i) = αn then |Λ 2 [M ]| = Λ 1 [M ] mαβ n mα 2 β, and
|E 1 ( Ẽ)| m(αβ) 2 1 - c 2 1 n E 2 ( Ẽ) m(αβ) 2 1 - c 2 1 n (2.16) Since m(αβ) 2 < 1 then 1 - c 2 1 n αβ-1
αβ(1-mαβ) . Finally we have

|E 1 ( Ẽ)| m(αβ) 2 αβ -1 αβ(1 -mαβ) , E 2 ( Ẽ) m(αβ) 2 αβ -1 αβ(1 -mαβ) ,
and then

|E 1 ( Ẽ)| αβ -1 1 -mαβ mαβ, E 2 ( Ẽ) αβ(αβ -1) 1 -mαβ mαβ.
(2.17)

The proof is now concluded.

It is worth stressing that our mathematical analysis has been limited to the case m(αβ) 2 < 1 because our numerical analysis has shown that the eigenvector V 1 [M ] and the degree vector M e are not likely to be close to each other when m(αβ) 2 1 .

An upper bound for m depending on the matrix elements

This section is concerned with the derivation of an upper bound for m. The advantage is that the upper bound of the error vectors (2.8) and (2.14) depends only on the elements of M . In order to obtain the main result of this section, we need some preliminary results. Proposition 2.2.1 Let M 1 , M 2 ∈ M n (R) be two symmetric matrices of size n. Then

λ i [M 1 ] + λ n [M 2 ] λ i [M 1 + M 2 ] λ i [M 1 ] + λ 1 [M 2 ], ∀i ∈ {1, 2, . . . , n} (2.18)
Proof The proof is based on paper [START_REF] Weyl | Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen[END_REF]. Since M 2 is a real symmetric matrix, then there exists an orthonormal base

B M 2 = {W 1 , W 2 , . . . , W n } of eigenvectors of M 2 whose spectrum is {λ 1 [M 2 ], λ 2 [M 2 ], . . . , λ n [M 2 ]}. If v ∈ R n , then v = n j=1 cj W j = (c 1 , c2 , ..., cn ) B M 2 .
Therefore

M 2 v = n j=1 cj λ j [M 2 ]W j , and
< M 2 v, v >= n j=1 (c j ) 2 λ j [M 2 ]. Since λ j [M 2 ] λ 1 [M 2 ]
, ∀j ∈ {1, 2, . . . , n}, we have

< M 2 v, v > λ 1 [M 2 ] n j=1 (c j ) 2 = λ 1 [M 2 ]||v|| 2 2 . (2.19)
We can assume that ||v|| 2 2 = 1 without loss of generality. We have

< (M 1 + M 2 )v, v >=< M 1 v, v > + < M 2 v, v > (2.19) < M 1 v, v > +λ 1 [M 2 ].
Let i ∈ {1, 2, . . . , n} and S ⊂ R n . Then

max dim(S)=i min x =0 x∈S < (M 1 + M 2 )v, v > max dim(S)=i min x =0 x∈S (< M 1 v, v >) + λ 1 [M 2 ].
Since M 1 and M 1 + M 2 are symmetric matrices, by the Min-Max theorem 1.1.3, we have

λ i [M 1 + M 2 ] λ i [M 1 ] + λ 1 [M 2 ].
By applying the previous inequality to

M 1 = M 1 + M 2 -M 2 , we have λ i [M 1 ] = λ i [M 1 + M 2 -M 2 ] λ i [M 1 + M 2 ] + λ 1 [-M 2 ] = λ i [M 1 + M 2 ] -λ n [M 2 ],
therefore the proof.

Remark 2.2.1 Let d i (M 1 , M 2 ) = λ i [M 1 + M 2 ] -λ i [M 1 ]. The proposition 2.2.1 states that λ n [M 2 ] d i (M 1 , M 2 ) λ 1 [M 2 ]. ( 2
(λ) = det(M -λ i [M ]Id -λId),
which admits λ = 0 as solution because λ i [M ] is an eigenvalue of M . Accordingly, the spectrum of Mi can decomposed as follows:

Spec( Mi ) = Spec + ( Mi ) Spec -( Mi ) Spec 0 ( Mi ),
where

Spec + ( Mi ) ={λ > 0, χ Mi (λ) = 0}, Spec -( Mi ) ={λ < 0, χ Mi (λ) = 0}, Spec 0 ( Mi ) ={λ = 0, χ Mi (λ) = 0} = ∅.
We assume that the cardinality of the sets Spec + ( Mi ), Spec -( Mi ), and Spec 0 ( Mi ) is s, t, and n -s -t, respectively. In particular ∃ k ∈ {1, 2, . . . , n} such that

λ 1 [ Mi ] • • • λ k[ Mi ] = 0 • • • λ n [ Mi ].
We can assume that k = i. Indeed, if V is the matrix whose columns are

V = (V 1 [M ] | V 2 [M ] | . . . | V n [M ]), then V -1 Mi V =V -1 (M -λ i [M ]Id)V =V -1 M V -V -1 λ i [M ]V =         λ 1 [M ] 0 . . . λ i [M ] . . . 0 λ n [M ]         -         λ i [M ] 0 . . . λ i [M ] . . . 0 λ i [M ]         =         λ 1 [M ] -λ i [M ] 0 . . . 0 . . . 0 λ n [M ] -λ i [M ]         . Since λ j [M ] -λ i [M ] 0, ∀j ∈ {1, 2, . . . , i -1} and λ j [M ] -λ i [M ] 0, ∀j ∈ {i + 1, . . . , n}, hence k = i.
According to the Sylvester theorem 1.1.2, D -1/2

Mi D -1/2 has s strictly positive eigenvalues, t strictly negative eigenvalues, and n -s -t zero eigenvalues. In particular, we have:

λ i [D -1/2
Mi D -1/2 ] = 0.

Setting D = D -1/2 Mi D -1/2 , we have

D = D -1/2 Mi D -1/2 = D -1/2 (M -λ i [M ]Id)D -1/2 = D -1/2 M D -1/2 -λ i [M ]D -1 , D (1.14) = D 1/2 M D -1/2 -λ i [M ]D -1 .

By proposition 2.2.1 with

M 1 = D 1/2 M D -1/2 , M 2 = -λ i [M ]D -1 , we have λ i [M 1 + M 2 ] = λ i [ D] = 0.
Then from Eq (2.20) we have

λ n [-λ i [M ]D -1 ] -λ i [D 1/2 M D -1/2 ] λ 1 [-λ i [M ]D -1 ]. (2.22) • Suppose λ i [M ] > 0. Then λ 1 [-λ i [M ]D -1 ] = -λ i [M ]λ n [D -1 ] = - λ i [M ] max j∈{1,2,...,n} (M e) (j) = - λ i (M )
αn . We also have

λ n [-λ i [M ]D -1 ] = -β λ i [M ]
n . Hence, by (2.22), we have

-β λ i [M ] n -λ i [D 1/2 M D -1/2 ] - λ i [M ] αn . Since λ i [D 1/2 M D -1/2 ] = λ i [M ],
and

-β λ i [M ] n -λ i [M ] - λ i [M ] αn < 0, we have |λ i [M ]| αn |λ i [M ]| β |λ i [M ]| n .
(2.23) (2.22) for the first equality), one has

• Suppose λ i [M ] = 0. Then λ i [D 1/2 M D -1/2 ] = λ i [M ] = 0 (from
λ i [M ] αn = λ i [M ] = β λ i [M ] n = 0.
(2.24)

• Suppose λ i [M ] < 0. Then λ 1 [-λ i [M ]D -1 ] = -β λ i [M ] n . We also have λ n [-λ i [M ]D -1 ] = - λ i [M ] nα . Then from (2.22) we have 0 < - λ i [M ] αn -λ i (M ) -β λ i [M ] n , with λ i [M ] < 0. Hence |λ i [M ]| αn |λ i [M ]| β |λ i [M ]| n . (2.25)
Then the proof is conluded.

Lemma 2.2.2 Let M ∈ SIR n (R + ), α and β defined by Eq (1.10) and Eq (1.11), and m defined by Eq (1.12). Let M = D -1 M be the matrix defined in (1.14). Then Proof Let N = ( t M ) be the transpose matrix of M . Then N and M have the same eigenvalues.

|λ i [M ]| m, ∀i ∈ {2,
In particular by Eq (1.16) we have

n i=1 N ij = n j=1 M ij = 1, ∀i ∈ {1, 2, . . . , n}.
It is easy to prove that the number 1 is an eigenvalue of M whose eigenvector is e (in particular t eN = t e). Moreover, the number 1 is the highest eigenvalue.

Let

B N = {V 1 [N ], V 2 [N ], . . . , V n [N ]
} be the set of the eigenvectors of N related to the eigenvalues

{λ 1 [N ], λ 2 [N ], . . . , λ n [N ]}. Let V (r) i [N ] = max j∈{1,2,...,n} |V (j) i |, for i ∈ {1, 2, . . . , n}, then |λ i [N ]V (r) i [N ]| = n j=1 N rj V (j) i [N ] n j=1 N rj |V (r) i | = |V (r) i |, hence λ i [N ] 1.
According to proposition 2 the matrix M is irreducible. Since M is non-negative, the Perron-Frobenius theorem implies that λ 1 R can be a null vector. We have

[N ] = 1 > |λ i [N ]|, ∀i = {2,
N V i [N ] = λ i [N ]V i [N ]. Hence N V i [N ] -0 = λ i [N ]V i [N ] ⇐⇒ (N -R t e)V i [N ] = λ i [N ]V i [N ]. It means that λ i [N ] is an eigenvalue of N -R t e with eigenvector V i [N ], for all i 2. Moreover (N -R t e) kj = N kj -R k = N kj -min j∈{1,2,...,n}
N kj 0, ∀k ∈ {1, 2, . . . , n}.

The Gerschgorin theorem 1.1.4 and the remark 1.1.4 imply that there exists j ∈ {1, 2, . . . , n} such that

|λ j [N ]| n k=1 |(N -R t e) kj | = n k=1 |N kj -R k | = n k=1 (N kj -R k ) n k=1 N kj - n k=1 R k 1 - n k=1 min l∈{1,2,...,n} M lk .
Let Q be the column vector whose k component reads:

Q k = max j∈{1,2,...,n} N kj = max j∈{1,2,...,n} M jk > 0, k ∈ {1, 2, . . . , n}.
One has

(Q t e -N )V i [N ] = λ i [N ]V i [N ]. It means that λ i [N ] is the eigenvalue of Q t e -N for the eigenvector V i [N ]. Note that Q k N kq , ∀q ∈ {1, 2, .
. . , n}. The Gerschgorin theorem 1.1.4 and the remark 1.1.4 imply that there exists j ∈ {1, 2, . . . , n} such that

|λ j [N ]| n k=1 |(Q t e -N ) kj | = n k=1 |Q k -N kj | = n k=1 (Q k -N kj ) n k=1 Q k - n k=1 N kj n k=1 max l∈{1,2,...,n} M lk -1.
Theorem 2.2.1 Let M ∈ SIR n (R + ), α and β defined by Eq (1.10) and Eq (1.11), and m defined by Eq (1.12). Let M = D -1 M be the matrix defined in (1.14)

k 0 ∈ N such that M k ∈ M n (R + ), ∀k k 0 .
We denote by EN n (R) the set of eventually nonnegative matrices.

It is straightforward that SIR n (R + ) ⊂ EN n (R). Theorem 2.3.1 ([48]) Let M ∈ EN n (R).
Then M possesses the Perron-Frobenius property, namely:

• Λ 1 [M ] > 0; • V (q) 1 [M ] 0, ∀q ∈ {1, 2, . . . , n}.
Bearing all above in mind, our result (Eqs (2.9) and (2.15)) can be generalized as follows.

Theorem 2.3.2 Let M ∈ EN n (R) be a diagonalizable matrix, α, β and m defined by Eq (1.10), Eq (1.11), and Eq (1.12), respectively. If |αβ| > 1 and | m|(αβ) 2 < 1 then

|E 1 (E)| | m| √ n (|αβ| -1) |αβ| (1 -| m|(αβ) 2 )(1 -| mαβ|) E 2 (E) | m| √ n |αβ| |αβ| -1 1 -| m|(αβ) 2 (2.32) and |E 1 ( Ẽ)| |αβ| -1 1 -| mαβ| | mαβ| E 2 ( Ẽ) |αβ|(|αβ| -1) 1 -| mαβ| | mαβ| (2.33)
where E and Ẽ are defined by Eq (2.8) and Eq (2.14), respectively.

On not diagonalizable matrix

Let M ∈ SIR n (R + ). It is easy to prove that [49] lim k→+∞ M k e ρ(M ) k = c 1 V 1 [M ]. (2.34) Indeed M k e = n i=1 c i (Λ i [M ]) k V i [M ], ∀k ∈ N * .
then

M k e = (Λ 1 [M ]) k c 1 V 1 [M ] + n i=2 c i Λ i [M ] Λ 1 [M ] k V i [M ] . Since Λ 1 [M ] is simple and Λ 1 [M ] > 0 then Λ 1 [M ] > Λ i [M ], ∀i ∈ {2, 3, . . . , n}, therefore lim k→+∞ M k e ρ(M ) k = c 1 V 1 [M ].
The aim of this section is to generalize the above result by supposing that M is not a diagonalizable matrix, namely M ∈ M n (R + ). Our results follows.

Lemma 2.4.1 Let M ∈ M n (R + ). Assume that there exists p n distinct eigenvalues ordered as follows

Γ 1 [M ] > ... > Γ p [M ]
and let m i be the algebraic multiplicity of Γ i [M ], for i ∈ {1, 2, . . . , p}, such that p i=1 m i = n. Then for all k 1, the following formula holds true:

M k = p i=1 m i -1 j=0 C j k Γ i [M ] k-j (M -Γ i [M ]Id n ) j Π i [M ], (2.35)
where

C j k = k! j!(k-j)! and Π i [M ] is the spectral projector on the ith-eigenspace related to Γ i [M ], ∀i ∈ {1, 2, . . . , p}.
Proof As known, the projector Π i [M ] satisfies the following properties [START_REF] Beck | Objectif Agrégation[END_REF]:

p i=1 Π i = Id n and Π 2 i = Π i , ∀i ∈ {1, 2, . . . , p}.
Since

M = M Id n = M p i=1 Π i [M ],
one has

M = p i=1 [Γ i [M ]Π i [M ] + (M -Γ i [M ]Id n )Π i [M ]] = p i=1 [Γ i [M ]Id n + (M -Γ i [M ]Id n )] Π i [M ],
and consequently

M k = p i=1 (Γ i [M ]Id n + (M -Γ i [M ]Id n )) k Π i [M ].
By using the Newton binomial formula, we have

(Γ i [M ]Id + (M -Γ i [M ]Id)) k =   m i -1 j=0 C j k Γ i [M ] k-j (M -Γ i [M ]Id n ) j   Id,
where we have used the property that the matrix

M -Γ i [M ]Id n , for i ∈ {1, 2, ..., p}, is m i - nilpotent [49].
The main result of this section follows.

Theorem 2.4.1 Let M ∈ M n (R + ). Assume that there exists p n distinct eigenvalues ordered as follows

Γ 1 [M ] > ... > Γ p [M ],
and the algebraic multiplicity of

Γ 1 [M ] is m 1 = 1. Then lim k→+∞ M Γ 1 [M ] k = Π 1 [M ], (2.36)
where Π 1 [M ] is the spectral projector on the eigenspace related to

Γ 1 [M ].
Proof By using Stirling formula [START_REF] Beck | Objectif Agrégation[END_REF] the following limit holds true:

lim k→+∞ C j k Γ i [M ] Γ 1 [M ] k = 0, ∀j ∈ {0, 1, . . . , m i -1}, ∀i ∈ {2, 3, . . . , p}. (2.37)
From Eq (2.35), we have

M Γ 1 [M ] k = p i=1 m i -1 j=0 C j k Γ i [M ] Γ 1 [M ] k M Γ i [M ] -Id n j Π i [M ]. (2.38)
From Eq (2.37) and since m 1 = 1 the only non-zero term, when k → +∞, of the previous equation is the one for i = 1.

An algebraic equation for Λ 1 [M ]

This section deals with the derivation of an algebraic equation for Λ 1 [M ].

Theorem 2.5.1 Let M ∈ SIR n (R + ). Then the highest eigenvalue Λ 1 [M ] is solution of the following polynomial equation:

x k+2 + < M k E, E > x 2 - < M k+2 e, e > c 2 1 = 0, ∀k ∈ N * . (2.39)
Proof From Eq (2.4) one has immediately that

< M k+2 e, e >= c 2 1 Λ 1 [M ] k+2 + n i=2 c 2 i Λ i [M ] k+2 , ∀k ∈ N * . (2.40) Moreover it is straightforward that n i=2 c 2 i Λ i [M ] k+2 = (c 1 Λ 1 [M ]) 2 < M k E, E >, ∀k ∈ N * , (2.41)
and then the proof.

It is worth noting that for k = 1, Eq (2.39) reads:

(Λ 1 [M ]) 3 + < M E, E > (Λ 1 [M ]) 2 - < M 3 e, e > c 2 1 = 0. (2.42)
The last equation allows us to say that when < M E, E >< 0 and if

| < M E, E > | <M 3 e,e> c 2 1 then Λ 1 [M ] tends to -< M E, E >,
which has a meaningful sense if E is seen as a vector which components are error parameters.

The assortativity of

M when |Λ 2 [M ]| is close to Λ 1 [M ]
This section deals with the analysis of the assortativity A D defined by Eq (1.34) when

|Λ 2 [M ]| is close to Λ 1 [M ]
. Specifically we will analyse the numerator A defined by Eq (1.35) and we will set some conditions under which the principal eigenvector is close to the degree vector componentwise.

Theorem 2.6.1 Let M ∈ SIR n (R + ) and A(M ) defined by Eq (1.35). Assume that

(i) ∃ > 0: Λ 2 [M ] = Λ 1 [M ] -; (ii) Λ 3 [M ] = Λ 4 [M ] = • • • = Λ n [M ] = Λ with Λ ∈ R * . Then A(M ) = 0 if and only if c 1 = √ n and c 2 = c 3 = • • • = c n = 0, where c i , ∀i ∈ {1, 2, . . . , n}, is defined by Eq (2.6). Proof First we prove that c 1 = √ n and c 2 = c 3 = • • • = c n = 0 imply A(M ) = 0. According to Eqs (2.2) and (2.4), we have e = c 1 V 1 and M e = c 1 Λ 1 V 1 . Then A = (< M e, e >)(< M 3 e, e >) -(< M 2 e, e >) 2 = (c 1 Λ 1 ) 4 (c 1 × c 3 1 -(c 2 1 )
2 ) = 0. Vice versa. From Eq (2.4) and the assumptions (i) and (ii) one has

M e = c 1 Λ 1 V 1 + c 2 (Λ 1 -)V 2 + Λ n i=3 c i V i ,
and then (up to the first order in )

< M e, e > = c 2 1 Λ 1 + (n -c 2 1 )Λ + c 2 2 (Λ 1 -Λ) -c 2 2 , < M 2 e, e > = c 2 1 Λ 2 1 + (n -c 2 1 )Λ 2 + c 2 2 (Λ 2 1 -Λ 2 ) -2Λ 1 c 2 2 + o( ), < M 3 e, e > = c 2 1 Λ 3 1 + (n -c 2 1 )Λ 3 + c 2 2 (Λ 3 1 -Λ 3 ) -3Λ 2 1 c 2 2 + o( ).
(2.43)

Bearing all above in mind, the function A reads:

A = A 0 + A 1 + A 2 2 + o( 2 ), (2.44)
where

(1) A 0 = c 4 1 Λ 4 1 + (n -c 2 1 )c 2 1 Λ 1 Λ 3 + c 2 1 c 2 2 Λ 1 (Λ 3 1 -Λ 3 ) + (n -c 2 1 )c 2 1 Λ 3 1 Λ + (n -c 2 1 ) 2 Λ 4 + (n -c 2 1 )c 2 2 Λ(Λ 3 1 -Λ 3 ) + c 2 1 c 2 2 Λ 3 1 (Λ 1 -Λ) + (n -c 2 1 )Λ 3 (Λ 1 -Λ)c 2 2 + c 4 2 (Λ 1 -Λ)(Λ 3 1 -Λ 3 ) -c 4 1 Λ 4 1 -(n -c 2 1 )c 2 1 Λ 2 1 Λ 2 -c 2 1 c 2 2 Λ 2 1 (Λ 2 1 -Λ 2 ) -(n -c 2 1 )c 2 1 Λ 2 1 Λ 2 -(n -c 2 1 ) 2 Λ 4 -(n -c 2 1 )c 2 2 Λ 2 (Λ 2 1 -Λ 2 ) -c 2 1 c 2 2 Λ 2 1 (Λ 2 1 -Λ 2 ) -(n -c 2 1 )c 2 2 Λ 2 (Λ 2 1 -Λ 2 ) -c 4 2 (Λ 2 1 -Λ 2 ) 2 = (n -c 2 1 ) c 2 1 (Λ 1 Λ 3 + Λ 3 1 Λ -2Λ 2 1 Λ 2 ) + c 2 2 (ΛΛ 3 1 -Λ 4 + Λ 3 Λ 1 -Λ 4 -2Λ 2 Λ 2 1 + 2Λ 4 ) + c 2 1 c 2 2 Λ 4 1 -Λ 1 Λ 3 + Λ 4 1 -Λ 3 1 Λ -2Λ 4 1 + 2Λ 2 1 Λ 2 + c 4 2 Λ 4 1 -Λ 1 Λ 3 -ΛΛ 3 1 + Λ 4 -Λ 4 1 -Λ 4 + 2Λ 2 1 Λ 2 = (n -c 2 1 )Λ 1 Λ(c 2 1 + c 2 2 )(Λ 1 -Λ) 2 -c 2 1 c 2 2 Λ 1 Λ(Λ 1 -Λ) 2 -c 4 2 (Λ 1 -Λ) 2 Λ 1 Λ = Λ 1 Λ(Λ 1 -Λ) 2 (n -c 2 1 )(c 2 1 + c 2 2 ) -c 2 1 c 2 2 -c 4 2 , = Λ 1 Λ(Λ 1 -Λ) 2 (c 2 1 + c 2 2 ) n i=3 c 2 i . Since 0 = Λ 1 > Λ = 0 and c 2 1 > c 1 > 1 > 0, then A 0 = 0 ⇐⇒ n = c 2 1 + c 2 2 . (2.45) (2) A 1 = -c 2 2 c 2 1 Λ 3 1 + (n -c 2 1 )Λ 3 + c 2 2 (Λ 3 -Λ 3 ) -3Λ 2 1 c 2 2 c 2 1 Λ 1 + (n -c 2 1 )Λ + c 2 2 (Λ 1 -Λ) + 4Λ 1 c 2 2 c 2 1 Λ 2 1 + (n -c 2 1 )Λ 2 + c 2 2 (Λ 2 1 -Λ 2 ) = -c 2 1 c 2 2 Λ 3 1 -(n -c 2 1 )c 2 2 Λ 3 -c 4 2 (Λ 3 1 -Λ 3 ) -3c 2 1 c 2 2 Λ 3 1 -3(n -c 2 1 )Λ 2 1 Λc 2 2 -3c 4 2 Λ 5 1 2(Λ 1 -Λ) + 4Λ 3 1 c 2 1 c 2 2 + 4Λ 1 Λ 2 c 2 2 (n -c 2 1 ) + 4Λ 1 c 4 2 (Λ 2 1 -Λ 2 ) = c 2 2 Λ(4Λ 1 Λ -Λ 2 -3Λ 2 1 )(n -c 2 1 -c 2 2 ) Since 4Λ 1 Λ -Λ 2 -3Λ 2 1 = 0 then we have A 1 = 0 ⇐⇒ n = c 2 1 + c 2 2 . (2.46) (3) A 2 = 3(n -c 2 1 ) 2 Λ 2 1 + 3(n -c 2 1 )nΛ 2 1 -4(n -c 2 1 ) 2 Λ 2 1 + 2n(n -c 2 1 )Λ 2 1 -(n -c 2 1 ) 2 Λ 2 1 + (n -c 2 1 )nΛ 2 1 = c 2 1 Λ 2 1 (n -c 2 1
), where we have used

n = c 2 1 + c 2 2 and t eM k e =c 2 1 Λ k 1 + (n -c 2 1 )Λ k 1 1 - Λ 1 k , ∀k ∈ {1, 2, 3}, t eM k e =nΛ k 1 -k(n -c 1 )Λ k-1 1 + k(k -1) 2 (n -c 1 )Λ k-2 1 2 + o( 2 ), ∀k ∈ {1, 2, 3},
and

< M e, e > = nΛ 1 -(n -c 2 1 ) , < M 2 e, e > = nΛ 2 1 -2(n -c 2 1 )Λ 1 + (n -c 2 1 ) 2 , < M 3 e, e > = nΛ 3 1 -3(n -c 2 1 )Λ 2 1 + 3(n -c 2 1 )Λ 1 2 + o( 2 ).
(2.47)

Then immediately A 2 = 0 ⇐⇒ n = c 2 1 . (2.48)
Bearing all above in mind we have

c 2 = c 3 = • • • = c n = 0. Remark 2.6.1 The condition c 2 = c 3 = • • • = c n = 0 implies that E = 0 R n , which means that (see Eqs (2.7) and (1.38)) V 1 [M ] = M e √
< M e, M e > .

Remark 2.6.2 If the hypothesis (i) and (ii) are replaced by the following one:

Λ 1 [M ] > Λ 2 [M ] = Λ 3 [M ] = Λ 4 [M ] = • • • = Λ n [M ] = Λ with Λ ∈ R.
Moreover, some straightforward numerical simulations allow to see that the theorem remains true if

(i) ∃ > 0: Λ 2 [M ] = Λ 1 [M ] -; (ii') Λ 3 [M ] ∼ Λ 4 [M ] ∼ • • • ∼ Λ n [M ],
where '∼' means that the eigenvalues are not too far from each other, or if (i) and (ii') are replaced by the following one:

Λ 1 [M ] Λ 2 [M ] ∼ Λ 3 [M ] ∼ Λ 4 [M ] ∼ • • • ∼ Λ n [M ]
, where ' ' means 'is much larger than' (the higher it is, the more the theorem is satisfied).

It is worth stressing that Eq (2.42) allows to obtain another proof of the fact that A(M ) = 0 if the error E = 0 R n (see Eq (2.8)). Indeed, from Eq (1.35) and Eq (2.42) one obtains

A(M ) = (Λ 1 [M ]) 3 c 2 1 (< M e, e >) -(< M 2 e, e >) 2 + (c 1 Λ 1 [M ]) 2 (< M e, e >)( t EM E). (2.49) If E = 0 R n Eq (2.4) implies that Λ 1 [M ] = < M e, e > c 2 1
and then

A(M ) = (Λ 1 [M ]) 3 c 2 1 (< M e, e >) -(< M 2 e, e >) 2 = c 2 1 Λ 3 1 (< M e, e > -c 2 1 M ) = 0. (2.50)

Applications to known graphs

This section is devoted to numerical simulations with the aim to investigate the role of αβ in the distance between the degree vector and the principal eigenvector. Specifically we will consider four different types of graph (see the first column of Figure 2.1) and their related adjacency matrices (second column), where we also depict the empirical histogram of the degree vector (third column), and we plot the degree vector versus the principal eigenvector (forth column).

In particular we will consider the following four graph types:

• Erdos-Reyni (ER) graph [START_REF] Erdos | On random graphs[END_REF]. This graph is obtained as follow. Starting with n unconnected distinct nodes, each pair of nodes is connected with probability p ∈ [0, 1]. When n is high, the average degree fluctuates around p, with standard deviation p(1 -p)/n. According to the central limit theorem, the histogram of the degrees approaches to a Gaussian. Thus when n increases the fluctuations of the components of the degree vector decreases, and then the distance between the degree vector and the principal eigenvector decreases.

• Scale-free (SF) graph [START_REF] Albert | Statistical mechanics of complex networks[END_REF]. The distribution of the degrees of the nodes follows a power-law γ. The realization of a SF graph is made by the barabasi game: with an initial connected network of n 0 n nodes, n -n 0 nodes are added one at a time, and each of these new nodes are connected to n 00 n 0 initial node(s) with a probability proportional to the initial degree of the n 00 nodes. The choices of n 0 and n 00 allow to control γ (typically γ 1). Differently from the ER graph, the histogram of the degrees is now not a Gaussian. The difference between the degree vector and the principal eigenvector is higher than that of ER's.

• Checkered matrix, or checkerboard (CB) graph. It is a graph splitted into two families.

Each node is mainly connected with nodes of the same family with a probability p as in the ER graph's case, but here with rare connexions with the nodes of the other family. Specifically, both families are independently set up, and then the whole network is rewired 100 times preserving the degrees. The whole graph is then connected. The size of both families is controled by a parameter k: the number of nodes in one family is n/2 + k (in the simulations n is even) and n/2 -k in the other family.

• Bipartite (BP) graph [START_REF] Asratian | Bipartite Graphs and their Applications[END_REF] graph. It is a graph also splitted into two families. Each node of each family is exclusively linked with all the members of the other family. Then we rewire 100 times the whole network preserving the degrees. We can then control the size of both families by varying k as for the CB graph's.

Specifically in this section, we plot the number

||R|| = nE 2 ( Ẽ) = n j=1 ( Ẽ(j) ) 2 = n j=1 1 (M e) (j) 2 n i=2 c i Λ i [M ]V (j) i [M ] 2 .
which is another error. The latter is a composite function of the error E 2 ( Ẽ) that we have considered in this chapter. In particular we consider now the error ||R|| because it is more sensitive to the size n of the matrix M . For the CB and BP graphs, perfect networks are represented by the highest and lowest dots respectively. For each of these dots, the network has been 'rewired' preserving the degrees 100 times: one degree-preserving rewire as the following procedure. Randomly take two linked pairs of nodes, say (a,b) and (c,d). Break the links, and randomly associate a to c or d. Then associate b to d or c. This has the effect to decrease the absolute value of the assortativity, vertically moving on the up right plot of figure 2.2. However, the assortativity is not the only graph parameter which is related to αβ and ||R||.

Looking at the CB graph, we see that the ability of the graph to be divided into modules changes the value of ||R||. When then to focus on the modularity. The modularity Q(G) of a graph G is the strength of a division of G into modules, and defined by the following equation:

Q(G) = n i,j=1 M ij < M e, e > - (M e) (i) (M e) (j) (< M e, e >) 2 δ(i, j), (2.51)
where δ(i, j) = 1 if nodes i and j belongs to the same module. The optimal modules are found by the algorithm detailled on the papar [START_REF] Pons | Computing communities in large networks using random walks[END_REF]: a random walk into the optimal modules must stay in it sufficiently longuely. The left panels of Figure 2.3 show the evolution of ||R|| with the modularity, and, for the CB and BP graphs, perfect networks are represented by the highest and lowest dots respectively. The rewire process has also been applied. On the left and right panels, for ER and SF graphs, the rewiring process makes one move up to down in a vertical line.

Applications to biology and finance

Finding the 3D-structure of DNA in the nucleus of cell has become a challenge over the past ten years. HiC experimental matrices (see chapter 3 for further analysis) are contact maps, which are adjacency matrices such that the (i, j)-element of the matrix is 1 if there is a physical contact between the ith and jth genes (the DNA is oriented), 0 otherwise. The left panel of the first row of Figure 2.4 shows an irreducible HiC matrix, and the center panel plots the principal eigenvector versus the degree vector. As this figure shows, the degree vector of the genomic matrix is an almost linear function of the highest eigenvector. As shown in the caption of Figure 2.4, the calculation of αβ and m should have led one to predict this behavior.

To make predictions of stock prices, we calculate the log-returns from the time series representing them: the log-return, or return, at time t is the logarithm of the price at time t out of the price at time t -1. In order to study the variability of a portfolio, one way is to construct the autocorrelation matrix, whose (i, j)-element is the time scalar product of asset returns ith and jth for a certain time duration. If the assets are all chosen from the same stock index, they will necessarily be correlated: the autocorrelation matrix is eventually non-negative. The left panel of Figure 2.4 shows the matrix returns associated with the stock index CAC 40. The center panel of Figure 2.4 plots the autocorrelation matrix from the returns. Finally, the right panel of Figure 2.4 plots the degree vector versus the highest eigenvector, and as this figure shows, the degree vector of the autocorrelation matrix is an almost linear function of the highest eigenvector. The reader is referred to the caption of Figure 2.4, to see that the calculation of αβ and m should have led one to predict this behavior. In a financial viewpoint, this might not be surprising: the investors trade on the index (here the CAC40), and all the assets composing the CAC40 are correlated with the CAC40. Then the assets are positively correlated (we say that they have a positive beta), and of course the matrix is positive. The closeness of the degree vector and the principal eigenvector is a characteristic of correlation, so as to the explosion of the highest eigenvalue.

A critical analysis on the pertinent literature

This section deals with a review with a critical analysis on the existing results on upper bounds for the error vectors (2.8) and (2.14), and the assortativity (1.34). The main aim is to compare the results of this chapter with the pertinent literature.

Mieghem et al. results

Let M be an irreducible adjacency matrix. Starting from the following equation

Λ 1 [M ]V (j) 1 [M ] = n q=1 M jq V (q) 1 [M ], ∀j ∈ {1, 2, . . . , n},
and under the assumptions that

• there exists a node j 1 with high degree;

• all the neighbors of j 1 have a high degree;

the authors of [START_REF] Mieghem | Double orthogonality and the nature of networks[END_REF] have proven that the magnitude of V

(j 1 ) 1 [M ] is large. Since Var[V 1 [M ]] = E 2 (V 1 [M ]) -(E 1 (V 1 [M ])) 2 , E 2 (V 1 ) = < V 1 [M ], V 1 [M ] > n = 1 n and E 1 (V 1 [M ]) = 1 n -Var[V 1 [M ]],
it is easy to see that when the variance Var[V 1 [M ]] increases (namely the assortativity increases) then the mean

E 1 (V 1 [M ]
) decreases. Moreover the authors of [START_REF] Mieghem | Double orthogonality and the nature of networks[END_REF] have shown that

Λ 1 [M ] = < M 2 e, e > < M e, e > and V 1 [M ] = M e √ < M e, M e > , (2.52)
The numerical simulations performed in [START_REF] Mieghem | Double orthogonality and the nature of networks[END_REF] show that the difference between the degree vector M e and the principal eigenvector V 1 [M ] decreases when the assortativity decreases. This property seems to be related to the fact that αβ tends to 1. Moreover Eq (2.52) holds true if and only if (M e) (i) = Λ 1 [M ], ∀i = {1, 2, . . . , n} (this is another property satisfied by a regular matrix). Contrary to [START_REF] Mieghem | Double orthogonality and the nature of networks[END_REF], we have found that when the assortativity decreases towards negative values the difference between the degree vector and the principal eigenvector increases (see section 2.7).

Recht et al. results

Let M be an irreducible adjacency matrix and ||V 1 [M ]|| p be the following norm:

||V 1 [M ]|| p =      (V (1) 1
[M ]) p + ... + (V (n) 1 [M ]) p 1 p if 1 p < +∞ max i∈{1,2,...,n} (V (i) 1 [M ]) if p = +∞
In [START_REF] Papendieck | On maximal entries in the principal eigenvector of graphs[END_REF] the authors have derived the following upper bound for the components of V 1 [M ]:

V (j) 1 [M ] max i∈{1,2,...,n} (V (i) 1 [M ]) Λ 1 [M ] p-2 1 + Λ 1 [M ] p-2 1 p
, ∀j ∈ {1, 2, . . . , n}.

(2.53)

In particular, for p = 2 the inequality (2.53) reads:

max i∈{1,2,...,n} (V (i) 1 [M ]) 1 √ 2 .
From Eq (2.6), we deduce that c 1 n/ √ 2. From Eq (2.2) we obtain

n = n i=1 c 2 i , and then c 1 √ n. It is easy to see that n √ 2 < √ n
for n < 2. Therefore our upper bound, c 1 √ n, appears more suitable.

Abueida et al. results

Let M ∈ M n (R + ) be an irreducible matrix and

s(M m ) := min i,j∈{1,2,...,n} M m ij max i,j∈{1,2,...,n} M m ij , m ∈ N * , (2.54)
where M m ij is the (i, j)-element of the matrix M m , for all i, j ∈ {1, 2, . . . , n}. The main interest in s(M m ) is concerned with its limit when m goes to +∞ [START_REF] Abueida | Inverse spread limit of a nonnegative matrix[END_REF]. The number s(M m ) has been analyzed in the paper [START_REF] Abueida | Inverse spread limit of a nonnegative matrix[END_REF], when the matrix M is symmetric and under the assumption that the algebraic multiplicity of Λ 1 [M ] is 1. Specifically, the authors of the paper [START_REF] Abueida | Inverse spread limit of a nonnegative matrix[END_REF] have shown that

lim m→+∞ s(M m ) = min 1 i,j n V (i) 1 [M ] V (j) 1 [M ] 2 ∈ R + .
(2.55)

In particular, as shown in [START_REF] Abueida | Inverse spread limit of a nonnegative matrix[END_REF], the above limit is

0 if M is nilpotent. It is worth noting that s(M m ) is different from (αβ) -1 . Indeed, one has (αβ) -1 def = min i∈{1,2,...,n} (M e) (i) max i∈{1,2,...,n} (M e) (i) = min i∈{1,2,...,n}   n j=1 M ij   max i∈{1,2,...,n} n k=1 M ik = s(M ).
Bearing all above in mind and since the formula (2.9) and (2.15) only depend on αβ, we state that s(M m ) is not appropriate for measuring the distance between the principal eigenvector V 1 [M ] and the degree vector M e considering our upper bounds. Therefore we will not be interested in further investigations on s(M m ).

Conclusion

This chapter deals with the closeness of two vectors associated to a nonnegative, symmetric and irreducible matrix. The first vector is the sum of each row of the matrix, column by column. This vector is the degree vector of a graph if the matrix is the adjacency matrix. The second vector is the eigenvector associated with the highest eigenvalue of the matrix. The idea comes from the fact that the highest eigenvalue is sometimes far from the others, and this characteristic can be found in many experimental matrices -especially in biology and finance. This chapter is a attempt of a rigorous formulation of the error made by considering that both vectors are equal.

Thus we have seen that the upper bounds of the first and second order moments of the error vector (being the componentwise substraction or division of both vectors) can only be expressed by the elements themselves of the matrix. The calculation of these upper bounds can be enough to predict the closeness of both vectors (and then the value of the highest eigenvalue), and then diagonalizing -highly time consumming for large data matrices -could be avoided. If the matrix is only eventually nonnegative, which is a generalization of nonnegative, symmetric and irreducible matrices, the results are still satisfied. At the same time, a nondiagonalizable matrix out of its spectral radius to the power of k, tends, when k tends to infinity, to the matrix's spectral projection on its highest eigenspace, which is a known result for diagonalizable matrices.

In parallel with these inovative results, the closeness between the two vectors can be observed in the Erods-Reyni networks. However, the Scale-Free networks play the antagonist role: the error vector's components are high. The checkerboard and bipartite networks (and any other network) play intermediary roles, and for some of them, the both vectors are close. As there are networks verifying or not the studied mathematical properties, distinct fields can be considered. On the one hand, in biology, the HiC matrix is nonnegative, symmetric and irreducible. The result of the present chapter is essential for the following one. On the other hand, in finance, the autocorrelation matrices are eventually nonnegative, and the developed theorems might lead to active portfolio management strategies.

Relating this property to more general matrices, even non-square, constitutes further research perspectives.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • •• • • • • • • • • • • • • • • • • • • •
0.02 0.04 0.06 0.08 0.10 0.02 0.04 0.06 0.08 0.10

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • •• 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Bipartite • • • • • • •• • • • • • • •• • • • • • • • •• • •• • •• • • • • • • • • • ••• • • • • • • • • • • • • • • • • •• •• •• • • •• • • • • • • • • • • • • • • • • •• • • • • • • ••• • • • • • • • • • • •• • • • • • • • • •• • • • • • • • • • • • •• • • • • • • • • • •• • • • • ••• • • •• • • • • • • •• • • • • • • • •• • • • • • • • • • •• •• • • • • • • • • •• •• •• • • • 0.04
0.06 0.08 0.10 0.12 0.04 0.06 0.08 0.10 0.12 q q q q q q q q q q q 0.5 2.0 3.5 0 n 0 initial node(s) with a probability proportional to the initial degree of the n 00 nodes. The 'concentration' of 1's decreases when we go down right in the adjacency matrix. The distribution of the degrees of the nodes follows a power-law γ (the loglog histogram is also shown), and the distance between the degree vector and the principal eigenvector is larger than that of an ER graph. Third row: the CB graph. It is a graph splitted into two families. Each node is mainly connected with nodes of the same family with a probability p as in the ER graph's case, with rare connexion with the nodes of the other family. We then rewire 100 times preserving the degrees. The adjacency matrix is a chessboard, and the degree distribution is as two Gaussians. The distance between the degree vector and the principal eigenvector increases when the difference of the size of both families increases. Fourth row: The BP graph. It is a graph also splitted into two families. Each node of each family is exclusively linked with all the members of the other family. We then rewire 100 times preserving the degrees. The histogram is far from being Gaussian and the error is large. Chapter 3

HiC Data Analysis

This chapter deals with the statistical analysis of HiC (High Chromosome contacts) matrices, which are genomic contact maps: the element M ij of the matrix is strictly positive if the ith and jth loci are in physical contact in the nucleus. The aim of the statistical study is to find the nuclear compartments from the analysis of a HiC matrix M . In section 3.1 we introduce the HiC concept. Section 3.2 gives a sketch of ideas why the highest eigenvalue is high for HiC matrices. Section 3.3 is the statistical analysis of the HiC matrix from human Embryonic Stem Cell (hESC) and also of the HiC matrix from human fibroblast. Finally section 3.4 deals with the 3D-reconstruction of DNA from a HiC matrix.

Introduction to HiC

With the development of DNA sequencing techniques, genomics has been trying to understand living beings by studying their genome. DNA in cell supports the genetic information, then characterizing each organism. In multicellular organism, cell differentiation is the set of processes which make a stem cell being specialized, e.g. fibroblast (connected tissue cell, for instance muscular). Then a natural interest arising is that the primary DNA structure (consecutive nucleotides sequence) cannot be responsible for differentiation by itself, because this structure is the same from the stem cell to the specialized cell. However the information about genome is different since the cell has changed: there are regions on DNA that are expressed, and others that are not anymore, and vice versa. But what makes a DNA region be expressed or not? An answer can be the proteins. Indeed, a group of proteins can be settled on a region of DNA to accelerate or to slow down the genetic expression. Moreover, proteins can bring different DNA parts in physical contact with each other, in order (i) to make the regions be expressed or not, and (ii) to strongly compact DNA in cell nucleus. In human being for instance, one DNA per nucleus of cell is indeed 2-meter long, whereas nucleus has a diameter of ∼ 1µm. All of this allows to argue that the 3D-structure of DNA in nucleus of cell depends on the specialization of the cell. In other words, the 3D-structure of DNA in stem cells must be different from the one in differentiated cells. Bearing all above in mind, the 3D-structure of DNA is linked to its biological function, yet our understanding of higher order genomic structure is rough and incomplete. In the nucleus of eukaryotic cells, chromosomes occupy distinct territories in the nucleus. Different regions from these chromosomes are also brought into physical contacts, thus forming nuclear compartments, and there exists numerous models that have been performed for showing how the chomosomes occupy the compartments. These models, however, provide only few mechanistic details about the relationship between higher order chromatin structure and genome function. In genomic technologies, some experimental revolutions in the study of 3D-structure of DNA have been employed. In particular, HiC has been introduced as a method for identifying chromatin interactions on the whole genome [START_REF] Lieberman-Aiden | Comprehensive mapping of long range interactions reveals folding principles of the humain genom[END_REF]. It has been proven that analyzing the HiC matrix can lead to the knowledge of nuclear compartments [START_REF] Lieberman-Aiden | Comprehensive mapping of long range interactions reveals folding principles of the humain genom[END_REF][START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF]. See Figure 3.2 for a simple example of the notion of nuclear compartment. For example, the two most known compartments for each cell is the one containing the euchromatin (not densed DNA) and the other one containing the heterochromatin (densed DNA) (chromatin = DNA + proteins), see Figure 3.2. It was largely admitted that euchromatin corresponds to expressed regions of DNA, whereas heterochromatin corresponds to not expressed regions of DNA. Here, there are two compartments: one corresponding to the biological feature 'expressed', and the other one 'not expressed'. One HiC experiment allows unbiased identification of chromatin interactions across an entire genome. Specifically cells are crosslinked with formaldehyde; DNA is digested with a restriction enzyme that leaves a 5'-overhang; the 5'-overhang is filled, including a biotinylated residue; and the resulting blunt-end fragments are ligated under dilute conditions that favor ligation events between the cross-linked DNA fragments. The resulting DNA sample contains ligation products consisting of fragments that were originally in close spatial proximity in the nucleus, marked with biotin at the junction. A library is created by shearing the DNA and selecting the biotincontaining fragments with streptavidin beads. The library is then analyzed using massively parallel DNA sequencing, producing a catalog of interacting fragments. From this catalog, a contact map M is produced, from which we can extract the information of the chromosome territories. The interested reader in a more deeper understanding of the HiC experimental procedure can refer to the paper [START_REF] Maze | Hic: a sequencing approach to characterize interactions over large genomic regions[END_REF]. Thus, in this chapter, a HiC matrix is considered to be a n × n symmetric and non-negative matrix M such that M ij = 1, if the ith DNA-segment and the jth DNA-segment are in physical contact, 0 otherwise. To find the compartments from the matrix M , a common and well-known technique is to bin the matrix M . For example, we usually reduce ten (resp. a hundred, etc.) consecutive segments into a unique segment which contains all the contact information from the initial segments: the order of the new matrix is n/10 (resp. n/100, etc.).

The literature proposes plenty of methods allowing to obtain the compartments. Accordingly, in paper [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF], the compartments are found in the intra-chromosomal interaction map, from the mouse ES cell genome. The process is based on a Hidden Markov Model: the observed directionality index characterizes the highly biased interaction frequencies at the periphery of the compartments. The boundary is then the place in the genome where the DI changes significativly. Each boundary is then identified. The Topological Associated Domains (TADs) are the domains in the chromosomes delimited by two consecutive boundaries. The paper [START_REF] Yaffe | Probabilistic modeling of Hi-C contact maps eliminates systematic biases to characterize global chromosomal architecture[END_REF] explained that there are three systematic biases in HiC experiments. The first is the restriction fragment length, correlated with contact probabilities. The second is the DNA GC content: it is shown that GC content is a source of incompatibilty between the replicates. The third is due to some interactions between non-specific cleavage sites. Taking into account these biases, a Bayesian model is performed to cluster a trans-contact map, knowing the biases. Then clustering the resulting matrix (by a k-means [START_REF] Hartigan | A K-means clustering algorithm[END_REF] procedure directly applied on the matrix) show three main columns, each of them corresponding to a genomic activity degree. Thus, they have identified (i) the boundaries, (ii) but also the compartments themself by their clustered matrix.

The paper [START_REF] Imakaev | Iterative correction of hi-c data reveals hallmarks of chromosome organization[END_REF] assumes that the biases are linear in the sense that the raw matrix element can be written as E ij = B i B j T ij . Here, B represents quantified biases (not necessarily the same as above, in [START_REF] Yaffe | Probabilistic modeling of Hi-C contact maps eliminates systematic biases to characterize global chromosomal architecture[END_REF]), calculated by maximum-likelihood estimations, and T is the raw matrix of contact map. We dispose of the biases by making a linear iterative correction using this formula. The final data matrix is smoother, do not have dashed lines, and is more homogene. After this correction, a Principal Component Analysis (PCA) [START_REF] Serre | Matrices: Theory and Applications[END_REF] is applied on the matrix: it consists of analysing the statistical scattering of the matrix elements over orthogonal axis. It is wellknown that the direction where the scattering is maximal is the direction of the first component, which is the direction of the eigenvector associated with the highest eigenvalue of the correlation matrix. The kmeans is used on the first and second eigenvectors, showing three distinct clusters.

The original first paper based on finding the compartments from HiC matrices is [START_REF] Lieberman-Aiden | Comprehensive mapping of long range interactions reveals folding principles of the humain genom[END_REF]. The authors calculated the empirical correlation matrix of the original HiC matrix: the element (i, j) of the resulting matrix is the empirical Pearson coefficient between the elements of the ith row and the elements of the jth column. This last matrix shows a very clear checkered pattern (exactly as the pattern of the adjacency matrix of a CB graph, see the previous chapter). This structure then suggests that each chromosome can be decomposed into two sets of loci (called A and B) such that contacts within each set are enriched and contacts between sets are depleted. The authors then partitioned each chromosome based on the sign of the eigenvector associated with the highest eigenvalue of the matrix. Specifically, the corresponding DNA parts having positive (negative) value of the eigenvector belong to the compartment A (B, respectively). They then compared the compartments to known genetic and epigenetic features. It is shown that the compartment A correlates strongly with the presence of genes and accessible chromatin.

In this chapter, we perform an analysis based on the two previous chapters of this thesis. Specifically a HiC matrix M , as an inter-chromosomal contact map, has itself the particularity to be a checkerboard (CB, see the previous chapter) matrix, see Figure 3.1. This property of HiC matrices allows to use clustering algorithms directly on M , not on its empirical covariance matrix. The eigenvectors show strong trends, and each trend corresponds to a boundary between two consecutive rectangles of the chessboard.

In the following, two problems will be addressed here:

• identifying compartments at high resolution from a contact map;

• deduce the 3D-structure of the chromosomes from this matrix.

A natural problem is enhanced from the first previous point: we do not know the number of compartments, which is not necessary 2. To find this number and bearing the chapters 1 and 2 in mind, we will not use the principal eigenvector (eigenvector associated with the highest eigenvalue of M ). Indeed, this one is redundant in the analysis, since it does not have any other information that the matrix itself. However, the other eigenvectors are relevant, in particular the second and the third ones. Thus, the idea is to do a clustering analysis in the 2D-space of these two eigenvectors in order to identify the distinct compartments. The previous second point underlines the following question: how can one have the information of the original DNA 3D-structure from the contact map? To answer this question, we propose a two-step algorithm, which name is ShRec3D, and which accommodates noisy contact matrices. The particularity is that this algorithm is much faster than already existing ones. Let v 1 (resp. v 2 ) be the number of non-zero elements in any row of type A (resp. B), except the diagonal. v 1 and v 2 can be interpreted as the number of contact(s) that any part of the DNA in A (resp. B) has with itself in the same compartment. We immediately have

The spectrum of a checkerboard matrix

n = v 1 + v 2 + 2. ( 3.1)
('+2' we count the diagonal element not included in v 1 (resp. v 2 ) in A (resp. B)). To calculate the determinant, we use the fact that one does not change the determinant of a matrix by changing a line (resp. a column by symmetry) by making a linear combination of lines (resp. a column by symmetry). The euchromatin compartment is the territory in the nucleus such that the DNA interaction with itself is low, and the heterochromatin compartment is the territory in the nucleus such that the DNA interaction with itself is high.

© pixfocus.com One has det(M ) = a a a a a a a a • • • a a a a a a a a • • • a a a a a a a a • • • a a a a a a a a • • • b b b b b b b • • • b b b b b b b • • • b b b b b b b • • • b b b b b b b • • • b b b b b b b • • • a a a a a a a a • • • a a a a a a a a • • • b b b b b b b • • • b b b b b b b • • • . . . .

.

We add to the first column the sum of the others, of the same type A, we factorize by a + v 1 a, and subtract the A rows with the first A row. One obtains

det(M ) = 1 a a a a a a a • • • 0 a -a 0 0 0 0 a a • • • 0 0 a -a 0 0 0 0 0 • • • 0 0 0 a -a 0 0 0 0 • • • b b b b b b b • • • b b b b b b b • • • b b b b b b b • • • b b b b b b b • • • b b b b b b b • • • 0 0 0 0 a -a 0 0 0 • • • 0 0 0 0 0 a -a 0 0 • • • b b b b b b b • • • b b b b b b b • • • . . . . .
det(M ) = 1 a a a a a • • • 0 a -a 0 0 0 0 • • • 0 a a -a 0 0 0 • • • 0 0 0 a -a 0 0 • • • b -b b b b b • • • 0 b -b 0 0 0 • • • 0 0 b -b b b • • • 0 0 0 b -b b • • • 0 0 0 0 b -b • • • 0 0 0 0 a -a 0 • • • 0 0 0 0 0 a -a • • • . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . × (a + v 1 a)(b + v 2 b).
We must calculate the determinant of a triangular matrix, which is easier than before. We deduce that det

(M ) = (a + v 1 a)(b + v 2 b)(a -a) v 1 (b -b) v 2 . ( 3
1 0 1 0 1 0 1 1 0 1 1 0 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 Figure 3.4:
Example of simple checkerboarded matrix, which is a kind of HiC matrix. In this example, there are two nuclear compartments A and B. The probability of contact in the compartment A is higher than between the two distinct compartments. An ideal case of HiC matrix would be that there is no any 1 in the AB and BA regions. is worth stressing that M and P M t P have the same eigenvalues, and the eigenvectors are also the same, up to a permutation P of their components:

M = n i=1 Λ i [M ] t V i [M ]V i [M ] ⇐⇒ P M t P = n i=1 Λ i [M ] t (P V i [M ])(P V i [M ]). ( 3
= (v 1 + 1 -λ)(n -(v 1 + 1) -λ)λ n-2 Sp(M ) = {v 1 + 1, n -(v 1 + 1), 0} (3.4)
Thus we immediately see that Λ

1 (M ) = max(v 1 +1, n-(v 1 +1)), Λ 2 (M ) = min(v 1 +1, n-(v 1 +1
)), and Λ 3 (M ) = • • • = Λ n (M ) = 0 (see the chapter 2, theorem 2.6.1). If the matrix is irreducible, some straightforward simulations show that the spectrum is approximatly the one shown by Eq (3.4).

Compartments identification in human Embryonic Stem Cell (hESC)

Finding the number of compartments remains an open question, and we give here one possible answer. It is possible that the methods are not robust or do not give the same result in another cell type.

3.3.1 Why V 2 and V 3 are good candidates for the clustering study • The Figure 3.6 shows a checkerboard matrix, with some 0's in the AA and BB squares, and some 1's in the AB and BA squares. The eigenvalues are also perturbed, but

Λ 1 (M ) ∼ Λ 2 (M ) and Λ 3 (M ) ∼ • • • ∼ Λ n (M ) ∼ 0 (see theorem 2.6.1 of chapter 2). Then V 1 is near M e/ √
< M e, M e >, even if there are two distinct clouds of dots in the plot of V 1 versus M e. Actually, these two clouds are each linked to a compartment A or B. In this example V 1 and V 2 show unsignificant patterns.

• If K is low, the checkered pattern disappears. In the matrix in Figure 3.7 there are not enough 1's in all the squares: this is a typical HiC matrix. The highest eigenvalue is far from the others and the difference between V 1 and M e is again low. One surprise is the pattern of V 2 and V 3 : it is the same as for the experimental HiC matrix (actually the other eigenvectors do not have any trend and look noisy).

As a conclusion, V 2 and V 3 seem to be good candidates for the statistical analysis.

The optimal number of compartments is found by informational criteria, based on kmeans-clustering

Let M be the hESC HiC matrix we want to analyse, and suppose M is sparse enough so that the checker pattern is not seen at all. We use the spectral decomposition [?] of M :

M = n i=1 Λ i [M ] t V i [M ]V i [M ].
This equation can be rewritten as:

M = n i=1 M i , where M i = Λ i [M ] t V i [M ]V i [M ].
As already discussed before, V 1 is close to the sum of each row of M , which is also noisy. As already mentioned, we do not analyse it for finding the compartments. The idea is that distinct clouds of points might be observed in the (V 2 ,V 3 ) plane. Each cloud is associated to a compartment, and a kmeans allows to better separate the clouds. HiC simulations: the matrix is a perturbed checkerboard, with a very low number of 1's. The highest eigenvalue is far from the others, and V 1 is still close to M e. In this case, V 2 and also V 3 show interesting trends, and might be analyzed in deeper details to find the checker pattern.

The kmeans algorithm used here is the one of Hartagan and Wang [START_REF] Hartigan | A K-means clustering algorithm[END_REF]. It divides the n points in the 2D-(V 2 ,V 3 ) plane into k clusters, such that the within-cluster sum of squares is minimized.

The particularity of this algorithm is that no movement of a point from a cluster to another will reduce the within-cluster sum of squares. However, one must enter the number of cluster(s) k so that the algorithm runs, whereas this is the parameter we need to know. This found parameter k must be optimal in a certain sense.

Actually, the optimality is performed by employing informational criterii. The Aikake informational criterion (AIC), the Bayesian informational criterion (BIC) and the Hannan-Quenn criterion (HQ) [START_REF] Akaike | A new look at the statistical model identification[END_REF][START_REF] Schwarz | Estimating the dimension of a model[END_REF] are given by:

AIC(k) = ln (σ (k)) + 2 k n , BIC(k) = ln (σ (k)) + k ln (n) n , HQ(k) = ln (σ (k)) + 2k ln (ln (n)) n .
Here σ (k) is the realized residual variance associated with the within-cluster sum of squares.

A minimum is expected for the three quantities, and the value of k associated to the minimal values of the above three criterii is the desired optimal number of clusters.

An arising problem is that if M is too much sparse then the trends in the vectors V 2 and V 3 will remain invisible. Specifically the order of magnitude of the sample sparsity of M is 10 -5 : the number of 1's in the matrix M is one hundred thousand times lower than n 2 . Such lack of information makes the analysis impossible. In particular, it is highly likely that there is at least one row whose elements are all zero, making M be a reducible matrix.

The Enrichment method is a way to increase the sample sparsity of the matrix M by adding 1's. There are two ways to enrich the matrix M :

• Enrichment 1 : the genomic neighbors of the HiC neighbors of each DNA segment are HiC neighbors. Some connexions are added in the graph, so that 1's appear in the matrix M above, below, on the left and right of all the elements having the value of 1.

• Enrichment 2 : the HiC neighbours of the HiC neighbours of DNA segment are HiC neighbors. This enrichment corresponds to the transformation M → M + M 2 , which in Graph Theory usually gives the number of paths of length 1 and 2. The zero element M ij is now equalled to 1 if there is a possible connection of length 2 in the original graph, between the nodes i and j.

Simulations show that combinaisons of enrichments (in practice 121 is sufficient) allow to find trends in V 2 and V 3 . Finally, the output of the kmeans algorithm allows to obtain a vector whose ith element is 1, 2, . . . , k if the ith corresponding DNA-region is in compartment A, B,..., respectively, for all i ∈ {1, 2, . . . , n}.

Not only do we have the number k of compartments, but also the location of these compartments on the whole genome.

Figure (3.8) fully sums up all the approach, and Figure (3.9) validates it with three simulated compartments A, B and C, where the number of clouds in the (V 2 ,V 3 ) 2D space is three. It is worth noting that the algorithm is capable of finding the small A part between C and B, which means that the precision is sufficient for finding small sizes of subrectangles among big ones in the matrix.

Eigen decomposition AIC(red) ; BIC(black) ; HQ(blue)

M 2 = Λ 2 t V 2 V 2 , whereas it is not clear for M 1 since V 1 is close to M e
, which is noisy. We find the optimal number of clusters k (here k = 2) thanks to the criteria AIC, BIC and HQ, by doing successive kmeans in the (V 2 ,V 3 ) space. k is then the number of compartment(s), and we obtain the location of each DNA fragment in one compartment.

Results on large experimental HiC contact map, and comparisons with biological features

After having validated our computational approach on simulated data, we determine the actual compartments in human Embryonic Stem Cell (hESC) using published HiC datasets (the matrix M has size 800, 000 × 800, 000, one fragment is ∼ 3, 000 base pairs (bp)) obtained from Dixon et al. [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF]. The associated HiC matrix represents only inter -chromosomal interactions. The problem with this matrix is that it is highly reducible, since there are plenty of rows and columns with zero elements. This matrix then can be binned, by aggregating contacts over multiplier restriction fragments (giving for instance the matrix of the Figure (3.1)). The problem is that binning the matrix reduces the resolution as the size of the matrix decreases. To keep the highest resolution of the ∼ 3, 000 bp for one matrix element, one can actually reduce the matrix by directly removing these rows and columns, but also the rows and columns entirely full of 1's, which, bearing the previous chapters 1 and 2 in mind, may constitute a bias for the analysis. We find that the optimal number of compartments k varies from 4 to 6. The enrichments 121 allows to clarify the true number, which is actually 4 (see the Figure Then there are in total five compartments: on the one hand A, B, C and D, and on the other hand X inf , corresponding to the eliminated fragments before the analysis (the ones responsible for rows and columns with zero elements). We computed the proportion of contacts found between each compartments (see the Figure Fragments belonging to the compartment X inf are found at the peri-centromeric regions of each chromosomes, suggesting that X inf is a compartment in which several centromers are gathered. Moreover the fragments belonging to the compartments B and C are found around the pericentromeric regions and the fragments belonging to the compartments A and D are found in centromere distal and telomere proximal regions, and are only found on long chromosome arms. Looking at the Figure (3.10C), it seems that there are a colocalization of the compartments A and D on the one hand, and of the compartments B and C on the other hand. Following the whole genome from the extremity 5' to the other extremity 3' fragment by fragment, Figure 3.11 shows a matrix m whose element m ij is the proportion of paths from the compartment i to the compartment j, with i, j ∈ {A, B, C, D, X inf }. Specifically, the element m ij is the ratio of the number of paths from the compartments i to j, out of the number of paths from the compartments i to A, i to B, and so forth (the resulted matrix obtained by following the whole genome from the extremity 3' to the other extremity 5' is the transpose of m). We observe that there are main paths from B, C and D to A, and from A to B. There are more paths between B and D than between B and C, and C and D. These two results confirm the proximity of compartments B and D, however compartment A occupies more the whole nucleus than the others. Indeed, the proportion of the compartments A, B, C, D and X inf in the whole genome is 0.3318 (highest), 0.2341, 0.1734, 0.2202 and 0.0405 respectively. However, this analysis has the downside to consider also the very fast changes: if we made a zoom in around exclusively light blue regions, we would likely find also many dark blue rays, which have been taken into account in the analysis. This misleading method then gives a result which does not correspond to what we interpret, looking at Figure (3.10C). We should perform an algorithm ignoring the rays to lead to the more intuitive colocalization interpretation. Matrix m whose element m ij is the proportion of paths from the compartment i to the compartment j obtained by following the whole genome from the extremity 5' to the other extremity 3' fragment by fragment, with i, j ∈ {A, B, C, D, X inf }. Specifically, the element m ij is the ratio of the number of paths from the compartments i to j, out of the number of paths from the compartments i to A, i to B, and so forth.

of A, B, C and D in the genome (sum in a row = 1). We have found no significant proportion related to the compartment X inf , so that we do not represent it. All of these features are taken from the UCSC genome browser website, peak ChiP-seq UCSC server, H1-hESC line cell (hg19).

Figure (3.12A) shows the replication timing and the Lamina Associated Domains (LADs) distributions over the compartments. The lamina protein is known to be linked with the heterochromatin (i.e. dense chromatin) and with the nucleus surface. The heterochromatin is dense and the corresponding DNA regions are not express. It is found that the proportion of the LADs is higher in the compartments A and B than in the compartments C and D, suggesting that the compartments A and B are inactive compartments, whereas C and D are active compartments: the genome loci belonging to the compartments C and D contain expressing genes. Thus, this characteristic of the compartments C and D (A and B) are to compare with the one of the compartment A (B) obtained in [START_REF] Lieberman-Aiden | Comprehensive mapping of long range interactions reveals folding principles of the humain genom[END_REF].

The replication timing (early/late type) is a quantification of first to last replicated DNA regions. Accordingly when the replication timing increases then the corresponding DNA loci are the first in the genome to be replicated. Bearing this information in mind, we found that the replication timing is higher in the compartments C and D than in the compartments A and B: the first replicated sequences are in the active areas.

Figure (3.12B

) depicts the proportion of binding sites for specific Transcription Factors (TFs) in each compartment (in a row, the sum is 1). A TF is a protein specifically binding on a DNA sequence and is responsible for the DNA transcription. Therefore the binding sites correspond to expressed genes and must be found in the active genomic areas. It is confirmed here that the majority of the TFs' binding sites are located in the loci belonging to the compartments C and D.

Another biological feature is the histone, which is a protein responsible for the DNA organization into structural units, called nucleosomes. Histone modifications are known to affect chromosome function through at least two distinct mechanisms: (i) altering the electrostatic charge of the histone resulting in a structural change; (ii) these modifications are binding sites for protein recognition modules. These characteristics suggest that histone marks (marking histone modifications on the genome) are interesting to compare with the genomic location of the compartments. Figure (3.12C) maps the proportion of histone marks in each compartment (here again, in a row, the sum is 1). A histone is a protein responsible for the DNA organization into structural units, called nucleosomes. The specific histone marks shown in this picture are mainly present in the active compartments C and D.

At this stage, there are several open questions because of counter-examples appearing in these results

. There are TFs, including Nanog, which are in the inactive compartment A. The TF Nanog is responsible for the specialization of stem cell [START_REF] Mitsui | The homeoprotein nanog is required for maintenance of pluripotency in mouse epiblast and ES cells[END_REF], so it might appear natural to find it in the inactive area in hESC. However it is found in massive active areas in mESC (mouse) [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF]. There are two histone marks which stand out from the others. On the one hand, the histone H2AZ is located in promoter regions [START_REF] Chunyuan | H3.3/H2A.Z double variant-containing nucleosomes mark 'nucleosome-free regions' of active promoters and other regulatory regions[END_REF], however in the mouse the histone H2AZ is merely absent in undifferentiated cells [START_REF] Rangasamy | Pericentric heterochromatin becomes enriched with H2A.Z during early mammalian development[END_REF]. Here, we find this last result: the histone H2AZ is found in inactive chromatin. On the other hand, H3K9me3 is as well in active areas as in inactive areas, probably referring to specific roles at the embryonic level: especially in mESC, histones are necessary for silencing a set of genes [64]. Figure (3.13) depicts the proportion of chromatin marks, also called chromatincolors, in each compartment (in a row, the sum is 1). The promoters, enhancers and insulators are mostly located in the fragments belonging to the compartments C and D, validating that these fragments are active genes. Moreover, the heterochromatin mark is mainly in the compartment A, which means that this compartment contains densed chromatin. This is in agreement with the fact that the compartment A occupancy in the genome is higher than the others, which is a characteritics of heterochromatin. It is worth pointing out that the part of the genome in the compartments C and D may contain the euchromatine and the part of the genome in the compartments A and B may contain the heterochromatin. Bearing the results about the TF Nanog and the histones H2AZ and H3K9me3 in mind, these last features are therefore located in the heterochromatin in hESC.

In addition, the repetitive DNA segments are significantly contained in the compartment X inf , and they are mainly located in the centromeric regions of chromosomes (these regions are called satellite DNA, see Figure 3.10).

Compartments identification in human Fibroblast

In

3D Reconstruction of genome from HiC contact map

The Fast-Floyd algorithm

We have seen that from the HiC contact map, one can extract information about the type of contact/interaction. We can go further by reconstructing the 3D structure of DNA. Several methods already exist [START_REF] Varoquaux | A statistical approach for inferring the 3D structure of the genome[END_REF] [66] ..., but the computational time is prohibitive for large data input.

We propose a two-step algorithm allowing to get the 3D structure. It is straightforward to get all the distances between a set of N points given their 3D coordinates (Fig 3 .15 a, step 1). From this matrix, one can also infer easily the binary contact matrix given a contact threshold (Fig 3 .15 a, step 2). The goal is to infer the optimal 3D coordinates knowing only the contact matrix. The mathematical problem of reconstructing a spatial structure from the distances between its elements is solved by distance geometry [START_REF] Havel | The theory and practice of distance geometry[END_REF]. These methods involve the computation of the first three eigenvectors in intermediary matrix, the Gram matrix (Fig 3 .15 a, steps 4 and5).

An important step in Multi Dimensional Scaling (MDS)-based methods of chromosome reconstruction is therefore the derivation of a complete set of distances from a (possibly sparse) contact map (Fig 3 .15 a, step 6). We introduce a weighted graph whose nodes are the N loci detected in the experiment. The length of a link is determined as the inverse contact frequency between its end nodes. We then take for the distance between any two nodes the length of the shortest path relating them on the graph, computed using the Floyd-Warshall algorithm. Our method accommodates binary contact maps (for example, single-cell HiC data) [START_REF] Nagano | Single-cell Hi-C reveals cell-to-cell variability in chromosome structure[END_REF] by taking link lengths equal to 1 between contacting points, or else infinite (no link). Although it is approximate and gives distances in a dimensionless unit, this shortest-path metric assigns a sound distance (symmetric and satisfying the triangular inequality) to all pairs of points. It offers a way to achieve the preprocessing step common to all 3C-based techniques of converting observed contact frequencies into a complete set of distances, independently of the downstream reconstruction method.

This algorithm, which we call 'shortest-path reconstruction in 3D' (ShRec3D), combines this shortest-path distance with MDS to achieve chromosome reconstruction (Fig 3 .15 a, steps 4 and 6). We tested the efficiency of ShRec3D in a controlled in silico case. We generated a yeast genome 3D structure represented as N = 26, 538 beads (each corresponding to approximately three nucleosomes) linked by springs accounting for intrachromosomal DNA connectivity. The 16 yeast chromosomes were confined into a nucleus of radius 1. We compared, for data sets of various sizes, both the reconstruction accuracy and the speed of ShRec3D and two other methods, BACH [START_REF] Hu | Bayesian inference of spatial organizations of chromosomes[END_REF] and ChromSDE [START_REF] Burton | Inference of spatial organizations of chromosomes using semi-definite embedding approach and hi-c data[END_REF]. All gave satisfactory results in terms of accuracy (Fig 3 .16 a); however, on a personal computer, the run time for our script ranged from tens of seconds for small data sets (∼1,000 points) to 50 h for the largest one (26, 538 points), several orders of magnitude faster than other methods (Fig 3 .16 b). The limiting step for ShRec3D computation time is the Floyd-Warshall algorithm computing shortest paths on the contact map, whose worst-case performance scales as O(N 3 ). We also found that the accuracy of the MDS reconstruction applied directly to inverse-frequency distances was poor (Fig 3 .16 2), demonstrating the importance of using our shortest-path metric before MDS reconstruction.

We tested and compared ShRec3D to the above-mentioned alternative methods in conditions closer to those of real Hi-C experiments. The robustness of the ShRec3D was tested by modifying the original binarycontact map and generated a controlled amount of disorder by moving randomly a given fraction of positive entries, and it was insured for noise levels lower than 1% (maximal level in typical Hi-C experiments (Fig 3 .16 c), see Online Methods of [START_REF] Lesne | 3D genome reconstruction from chromosomal contacts[END_REF]). The probabilistic nature of BACH makes it efficient in the presence of high levels of noise; however, it remains slower than ShRec3D reconstruction by several orders of magnitude and is thus limited to small-size structures (Fig 3 .16 b). We then reproduced the superposition of single-cell contact maps reflecting the genome fold variations over a cell population that is characteristic of Hi-C maps; accordingly, one could reach only an average 3D structure. From a Langevin dynamic simulation [START_REF] Hajjoul | High-throughput chromatin motion tracking in living yeast reveals the flexibility of the fiber throughout the genome[END_REF] of our in silico genome, we extracted a variable number k of independent structures and computed the average of their contact maps (Online Methods of [START_REF] Lesne | 3D genome reconstruction from chromosomal contacts[END_REF]). The distances reconstructed with ShRec3D from this simulated average Hi-C contact map quantitatively matched the average distances in the superposition of structures (Fig 3 .16 d). This was also achieved by the alternative methods for a large number of structures; however, the comparison had to be limited to coarse-grained structures with 480 points, the largest size manageable in a reasonable time using BACH and ChromSDE. The increase in quality of MDS applied to inverse-frequency distances with the number of structures was expected, as the inverse-frequency expression becomes closer to the shortest-path distance when the average contact map becomes denser.

We implemented ShRec3D on experimental Hi-C data obtained in human embryonic stem cells [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF] and lymphoblastoids [START_REF] Kalhor | Genome architectures revealed by tethered chromosome conformation capture and population-based modeling[END_REF], exploiting both the very sparse Hi-C data obtained at the best available genomic resolution (restriction fragments) and coarse-grained data sets (where loci correspond to many restriction fragments). ShRec3D's ability to visualize average structures at different scales is illustrated by reconstructing a 30-Mbp region of chromosome 1 at 3-kb resolution (Fig 3 . d-f) would make it an efficient tool for genome scaffolding from Hi-C data [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] [START_REF] Kaplan | High-throughput genome scaffolding from in vivo dna interaction frequency[END_REF]. Alternative methods (BACH, MDS applied to inverse frequency distances and ChromSDE [START_REF] Zhang | 3d chromosome modeling with semi-definite programming and Hi-C data[END_REF]) did not manage to properly reconstruct fine-resolution structures in reasonable amounts of time. The potential of ShRec3D to devise 3D genome browsers is illustrated with the coloring of a 3D structure of chromosome 1 at resolution 30-kb according to the chromatin partition in two compartments. Any chemical, structural or functional annotation available on linear genomes can be similarly overlaid on chromosome 3D structures (for example, two histone H3 modifications, H3K9Ac and H3K9me3 (GEO GSM409308)). ShRec3D involves no ad hoc constraints or tunable parameters and is free from convergence issues and misleading transient outcomes. Its speed makes it applicable to both 3C or carbon-copy 3C (5C) data sets, which typically involve tens of loci, and high-resolution Hi-C data sets, comprising sparse contacts between hundreds of thousands of points. Its accurate reconstruction of average distances between genomic loci and visualization of a consensus structure enable a meaningful use of cell-population Hi-C data, especially when extended into 3D genome browsers.

The DNA 3D-structure of chromosome 1 and the compartments

In this conclusive subsection, The Fast-Floyd algorithm and the compartimentation are both taking into account for the chromomosome 1 of hESC. Specifically we apply the previous Fast-Floyd algorithm to the hESC intra-chromosomal HiC matrix, with the knowledge of the distribution of the DNA fragments in the nuclear compartments A, B, C and D obtained from the hESC inter -chromosomal HiC matrix already studied in this chapter. Figure 3.18 is a representation of the chromosome 1 in the hESC nucleus. The color code is the same of that of Figure 3.10C. The big red ball is the centromer (associated to the compartment X inf ), and the two black balls delimit the chromosome. We finally observe on this chromosome that the nuclear compartments in which are located the different parts of the chromosome 1 are globally well distinct. This is in accordance with the initial idea refered in paper [START_REF] Lieberman-Aiden | Comprehensive mapping of long range interactions reveals folding principles of the humain genom[END_REF] of nuclear compartimentation.

Conclusion

This chapter deals with data analysis of related chromosomal conformation for the whole humain Embryonic Stem Cell (hESC) genome. These data allows to construct frequency matrices at the resolution of 3k-base pairs, such that the element (i, j) of the matrix is zero if there is no any read of physical contact between the loci ith and jth of the genome. The interesting point of such matrices is that they are checkerboard matrices, studied in the previous chapters. The goal is then to identify nuclear compartments from this noisy matrix, and to give a biological sens to them. Differently from the previous litterature, we do not analyze the empirical correlated matrix, but the matrix itself. Because the matrix is too sparse, we removed the empty rows/columns, then giving one default compartment X inf . We found from the resulting irreducible matrix that the sum of each row, column by column, of the matrix is very close to the eigenvector associated to the highest eigenvalue. From the previous chapter, and since the matrix -hence the sum of the rows -is not very informative, the second and third eigenvectors are analyzed. These eigenvectors have the ability to have trends. Specifically, a kmeans is performed on the scatter plot of the two lattest eigenvectors, where distinct clouds of dots are revealled -a clear characteristic of checkerboard matrices. Informational criterii confirmed that there are 4 compartments, A, B, C and D, and 5 with X inf . The interactions are high between A and B and between C and D. Moreover, A and D are colocalized around the center of the chromosomes whereas B and C are colocalized mainly at their peripheries. The biological features considered in this thesis allowed to give a sense to these compartments. Thus, the replication timing and the Lamine Associated Domains (LADs) led to think that A and B are inactive compartments, whereas C and D are active compartments. These characteristics are confirmed by the histone marks and the Transcription Factors (TFs), which are mainly located in C and D. Some of the ways to confirm that C and D are active might be increased in efficiency. The compartments can give an idea about the 3D structure of DNA in the nucleus. Specifically, from a HiC matrix, we performed an algorithm, which allows to give the 3D structure of DNA in an efficiency running time (decreasing of two magnitude orders in comparison with existing algorithms). Thus, from the HiC matrix (seen as an adjacency matrix), we estimated the Euclidean distances between the loci by the Fast-Floyd algorithm. Some numerical simulations have shown that the recombined DNA 3D structure is very close to the original 3D structure. We have performed the analysis to the Fibroblast cell, and 2 compartments have been found. More biological features would lead to the identification of the compartments for each cellular type. These considerations would aim at a higher global understanding of cellular differentiation. Moreover, analyzing the HiC matrix associated to tumorous cells would also make us understand more about cancer.

Part II

Development and Analysis of Models in Biology

Chapter 4

A Delayed Mathematical Model for the Acute Inflammatory Response

This chapter deals with further developments on a mathematical model recently proposed for the modeling of the acute inflammatory response to infection or trauma. In particular in order to take into account that some interactions have not an immediate effect, we introduce time delays. Specifically the chapter deals with the existence of steady states, determining the parameter regimes where the equilibrium points are stable, and the onset of Hopf bifurcation appears. Numerical simulations are performed with the main aim of supporting the analytical results and investigate further dynamics. It is worth pointing out that the contents of this chapter constitutes the paper [START_REF] Bianca | A delayed mathematical model for the acute inflammatory response to infection, to appear[END_REF].

Introduction

The problem of abnormal organ repair has gained much attention considering that there is a significant shortage of organs available for transplantation. In this context the normal repair process, i.e. wound healing process, assumes an important role. Wound healing is an complex process by which the skin or organ repairs itself after injury [START_REF] Kirsner | The wound healing process[END_REF][START_REF] Velnar | The wound healing process: an overview of the cellular and molecular mechanism[END_REF]. Specifically wound healing comprises three sequential, overlapping, phases: the inflammation phase (hemostasis and the actual inflammation), the proliferation phase, and the maturation (remodeling) phase. Hemostasis occurs immediately after tissue injury and can be compared with the acute phase reaction of the innate immune system during infection. The first cells to appear in the wound area are neutrophils which start with the critical task of phagocytosis in order to destroy and remove bacteria, foreign particles and damaged tissue. Phagocytotic activity is crucial for the subsequent processes, because acute wounds that have a bacterial imbalance will not heal. The macrophage becomes the predominant inflammatory cell type in clean noninfected wounds. Every phase of the healing process consists of complex interactions between cells and mediators which tend to regulate the process. Cells participating in wound healing must be activated, i.e. undergo phenotypic alterations of cellular, biochemical, and functional properties.

During the inflamation phase, the immune system performs a fundamental action, see [START_REF] Cooper | Evolution of immune systems from self/not self to danger to artificial immune systems (AIS)[END_REF][START_REF] Parkin | An overview of the immune system[END_REF][START_REF] Alam | A brief review of the immune system[END_REF]. The response of the immune system to an infectious agent is subdivided into two main categories: Innate (non-specific) immunity response, which is mediated by granulocytes, macrophages, and NK cells [START_REF] Mackay | Innate immunity[END_REF]; Adaptive (specific, acquired) immune response, which is mediated by the lymphocytes [START_REF] Sallusto | Two subsets of memory t lymphocytes with distinct homing potentials and effector functions[END_REF]. The innate immune system is constitutively active and reacts immediately to infection. The adaptive immune response to an invading organism takes some time to develop.

Different mathematical models have been proposed for the modeling of immune system response [START_REF] Schlub | Comparing the kinetics of nk cells, CD4, and CD8 T cells in murine cytomegalovirus infection[END_REF][START_REF] Kumar | The dynamics of acute inflammation[END_REF]. Specifically mathematical models based on ordinary differential equations [START_REF] Nieman | A two-compartment mathematical model of endotoxin-induced inflammatory and physiologic alterations in swine[END_REF][START_REF] Paul | Modelling and analysis of time-lags in some basic patterns of cell proliferation[END_REF][START_REF] Pennisi | A mathematical model of immune-system-melanoma competition[END_REF][START_REF] Bianca | Mathematical modeling of the immune system recognition to mammary carcinoma antigen[END_REF], partial differential equations [START_REF] Eftimie | Interactions between the immune system and cancer: a brief review of non-spatial mathematical models[END_REF][START_REF] Waniewski | Mathematical Models of the Immune Response[END_REF], kinetic theory approach [START_REF] Kolev | Mathematical modelling of the competition between tumors and immune system considering the role of the antibodies[END_REF][START_REF] Bianca | A thermostatted kinetic framework with particle refuge for the modeling of tumor hiding[END_REF] and continuum mechanics approach [START_REF] Mollica | Mechanics in tumour growth[END_REF]. In the pertinent literature computational models have been also proposed, see [START_REF] Pappalardo | Induction of T-cell memory by a dendritic cell vaccine: a computational model[END_REF] and the review paper [START_REF] Pappalardo | Mathematical and computational models in tumor immunology[END_REF]. However the previous cited models are based on instantaneous interactions thus avoiding to take into account that various phenomena occur with some delay. In order to overcome this issue, delayed models have been proposed in population dynamics [START_REF] Cunningham | Time Lag in Population Models[END_REF][START_REF] Bianca | Modeling complex systems with particles refuge by thermostatted kinetic theory methods[END_REF], in immunology [START_REF] Marchuk | Mathematical Modelling of Immune Response in Infectious Diseases[END_REF][START_REF] Forys | Interleukin mathematical model of an immune system[END_REF], for tumor formation [START_REF] Byrne | The effect of time delays on the dynamics of avascular tumor growth[END_REF][START_REF] Piotrowska | Hopf bifurcation in a solid avascular tumour growth model with two discrete delays[END_REF][START_REF] Shi | Hopf bifurcation of a mathematical model for growth of tumors with an action of inhibitor and two time delays[END_REF], for economic systems [START_REF] Invernizzi | On lags and chaos in economic dynamic models[END_REF][START_REF] Zak | Kaleckian lags in general equilibrium[END_REF][START_REF] Bianca | The Cai model with time delay: existence of periodic solutions and asymptotic analysis[END_REF][START_REF] Gori | Hopf bifurcation in a cobweb model with discrete time delays[END_REF] and specifically in the study of the Solow model [START_REF] Bianca | Existence of limit cycles in the solow model with delayedlogistic population growth[END_REF][START_REF] Ferrara | Nonlinear dynamics in a solow model with delay and non-convex technology[END_REF], of the Dalgaard-Strulik model [START_REF] Bianca | On the dalgaard-strulik model with logistic population growth rate and delayed-carrying capacity[END_REF], of the credit risk contagion [START_REF] Bianca | Hopf bifurcations in a delayed microscopic model of credit risk contagion[END_REF], and in the asset of price [START_REF] Dibeh | Speculative dynamics in a time-delay model of asset prices[END_REF]. The introduction of time delay allows to enrich the description of the dynamics of a system in particular by changing the stability of the steady state and triggering onset of Hopf bifurcation.

This chapter is concerned with further investigations on a mathematical model recently proposed in [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] for the modeling of the acute inflammatory response to infection or trauma. In particular in order to take into account that some interactions have not an immediate effect, we introduce time delays. Specifically the chapter deals with the existence of steady states, determining the parameter regimes where the fixed points are stable, and the onset of Hopf bifurcation appears. As known a steady state belongs to the nullclines of the system and it is stable if the real part of each eigenvalue associated with the linearized system at that fixed point is negative. A bifurcation occurs when a change in a parameter alters the number of fixed points and/or their stability.

The contents of the present chapter are organized as follows. After this introduction, Section 6.2 is devoted to the mathematical analysis (existence of steady states and Hopf bifurcation) of the delayed pathogen equation, which shares various properties with the delayed model proposed in Section 3 which consists of a system of two delayed differential equations where the independent variables represent the levels of pathogen and the activated phagocytes (e.g. neutrophils); the model is built up from consideration of direct interactions of fundamental effectors and does not include components of the adaptive immune response, i.e. T-cells and specific antibodies. Numerical simulations are also performed within the Section 6.3 with the main aim of supporting the analytical results and investigating further nonlinear dynamics. Finally Section 6.4 is concerned with concluding remarks and further research perspectives.

The delayed pathogen equation

It is known that when an individual undergoes an injury, at the wound there exists a sentinel level of immune system cells (local immune response M ) able to respond and remove local infections (pathogens P ). According to [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] we assume that:

• M is inhibited at the rate k mp when it interacts with P at time t;

• P is inhibited at the rate k pm when it interacts with M at time t;

• s m models a source of M ;

• µ m models the death of M . Moreover we introduce a time delay τ ≥ 0 into the equation of P as follows:

     . M = s m -µ m M -k mp M P . P = -k pm M P τ
where P τ = P (t -τ ). Bearing [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] in mind thus we have

. P = k pg P 1 - P p ∞ - k pm s m P τ µ m + k mp P , (4.1)
where k pg is the pathogen growth rate and p ∞ is the carrying capacity of the pathogen population.

In what follows, the above parameters will be assumed non-negative.

Steady states and stability analysis

The steady states of Eq. (4.1) are such that the time derivative . P vanishes identically. It is immediate to see that the fixed points of Eq. (4.1) coincide with those for τ = 0. In particular P * = 0 is always a fixed point of Eq. (4.1). The other fixed points P * are solution of the following algebraic equation:

k pg 1 - P * p ∞ - k pm s m µ m + k mp P * = 0, namely if k pg k mp P 2 * + k pg (µ m -k mp p ∞ ) P * + p ∞ (k pm s m -k pg µ m ) = 0. (4.2)
The discriminant of Eq. (

= k pg k pg (µ m + k mp p ∞ ) 2 -4k mp k pm s m p ∞ , 4.2) reads ∆
then Eq. (4.1) also admits the fixed point

P * = (k mp p ∞ -µ m ) /2k mp if k pg = 4k mp k pm s m p ∞ /(p ∞ k mp + µ m ) 2 ,
and the fixed points P * = p 1 and P * = p 2 , p 1 < p 2 , with

P * = k pg (k mp p ∞ -µ m ) ± √ ∆ 2k pg k mp if k pg > 4k mp k pm s m p ∞ /(p ∞ k mp + µ m ) 2 .
In [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] Reynolds et al. have shown that P * = 0 is stable for k pg < k pm s m /µ m , and P * = p 2 is stable whenever it exists. Henceforth, in what follows we will deal only with these two stable fixed points.

Bearing all above in mind, the linearized equation of Eq. (4.1) around one of these two stable fixed points P * is . P = a (P -P * ) + b (P τ -P * ) ,

where if P * = 0 we have

a = k pg and b = - k pm s m µ m and if P * = k pg (k mp p ∞ -µ m ) + √ ∆ 2k pg k mp a = k pg 1 - 2P * p ∞ + 1 - P * p ∞ 2 k mp k pg P * k pm s m and b = -k pg 1 - P * p ∞
The associated characteristic equation of (4.4) reads λ -a -be -λτ = 0. (

It is well known that the fixed point P * of Eq. (4.1) is locally asymptotically stable if each of the characteristic roots of Eq. (4.5) has negative real parts. Hence, the marginal stability is determined by the equations λ = 0 and λ = iω, ω > 0. It is clear that the case λ = 0 cannot occur because a + b < 0. Let λ = iω be a root of the characteristic equation (4.5) with ω > 0. Substituting it into (4.5) and separating the real and imaginary parts yields

a = -b cos ωτ, ω = -b sin ωτ. (4.6)
Squaring each equation in (4.6), taking the sum and employing sin 2 ωτ + cos 2 ωτ = 1, we have

ω 2 = b 2 -a 2 . (4.7)
It is easy to see that Eq. ( 4.7) has one positive solution

ω 0 = b 2 -a 2
if |b| > |a|. From (4.6), one can obtain the value τ 0 corresponding to ω 0 as follows:

τ 0 =          1 ω 0 tan -1 ω 0 a , if a > 0, 1 ω 0 tan -1 ω 0 a + 2π, if a < 0.
One can also see that the purely imaginary root iω 0 is simple. If we suppose by contradiction λ = iω 0 to be a repeated root of (4.5), then differentiating (4.5) with respect to λ, inserting λ = iω 0 , and using (4.5), leads to ω 0 = 0, which gives a contradiction.

On the Hopf bifurcation

The conditions under which a Hopf bifurcation occurs at τ 0 are verified except for the transversality condition. Let λ (τ ) be the root of (4.5) near τ = τ 0 such that Re(λ(τ 0 )) = 0 and Im(λ(τ 0 )) = ω 0 . Differentiating both sides of Eq. (4.5) with respect to τ, we have dλ dτ

-1 = - 1 λ(λ -a) - τ λ .
Thus, we obtain

sign d (Reλ) dτ τ =τ 0 ω=ω 0 = sign Re dλ dτ -1 τ =τ 0 ω=ω 0 = sign 1 ω 2 0 + a 2 . (4.8)
Since the sign of (4.8) is positive, when λ = iω 0 , the only crossing of the imaginary axis is from left to right as τ increases. Consequently, the stability of the steady state P * can only be lost and not regained.

Bearing all the above analysis in mind, we can state the main result of this section.

Theorem 4.2.1 Let ∆ be defined by Eq. (4.3).

1) If k pg < k pm s m /µ m , then the steady state P * = 0 of Eq. (4.1) is locally asymptotically stable for τ < τ 0 and unstable for τ > τ 0 , where

τ 0 = 1 ω 0 tan -1 ω 0 k pg , ω 0 = k pm s m µ m 2 -k 2 pg .
Furthermore, Eq. (4.1) undergoes a Hopf bifurcation at P * = 0 when τ = τ 0 .

2)

If k pg > 4k mp k pm s m p ∞ /(p ∞ k mp + µ m ) 2 and
1 - P * p ∞ > 1 - 2P * p ∞ + 1 - P * p ∞ 2 k mp k pg P * k pm s m ,
the steady state

P * = p 2 = k pg (k mp p ∞ -µ m ) + √ ∆ /2k pg k mp
of Eq. (4.1) is locally asymptotically stable for τ < τ 0 and unstable for τ > τ 0 , where

τ 0 =          1 ω 0 tan -1 τ0 , if α > 0, 1 ω 0 tan -1 τ0 + 2π, if α < 0, and τ0 = ω 0 k pg 1 - 2P * p ∞ + 1 - P * p ∞ 2 k mp k pg P * k pm s m , α = 1 - 2P * p ∞ + 1 - P * p ∞ 2 k mp k pg P * k pm s m , ω 0 = k pg 1 - P * p ∞ 2 -α 2 .
Furthermore, Eq. (4.1) undergoes a Hopf bifurcation at P * = p 2 when τ = τ 0 .

The delayed model with the immune system response

In this section we couple the pathogen equation analyzed in the previous section with the role of phagocytic immune system cells (neutrophils and macrophages). Bearing paper [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] in mind, we modify the model proposed in [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] by introducing two time delays τ 1 ≥ 0 and τ 2 ≥ 0 as follows:

.

P = k pg P 1 - P p ∞ - k pm s m P τ 1 µ m + k mp P -k pn N * P, (4.9)
.

N * = s nr k nn N * τ 2 + k np P τ 2 µ nr + (k nn N * + k np P ) -µ n N * , (4.10)
where P and N * represent the levels of pathogen and activated phagocytes, respectively and • k pn is the rate at which N * consumes P ;

• s nr is the source of resting phagocytes;

• µ nr is the decay rate of resting phagocytes;

• k nn is the activation of resting phagocytes by previously activated phagocytes and their cytokines;

• µ n is the decay rate of activated phagocytes;

• k np is the activation rate of resting phagocytes by pa thogen.

It is worth stressing that with respect to paper [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF], the above model is derived by considering the delayed equations

. P = -k pm M P τ 1 and . N * = (k nn N * + k np P ) τ 2 N R -µ n N *
, where N R is the population of the resting phagocytes.

Steady states and stability analysis

The steady states of system (4.9)-(4.10) are obtained by setting . P = . N * = 0, P τ 1 = P and N * τ 2 = N * for all t. Therefore, when there is no time delay, i.e. τ 1 = τ 2 = 0, we recover the model considered in Section 2.2 of paper [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF], and in particular there exists the fixed point (P, N * ) = (0, 0) (health steady state). It is important to note that the delayed model admits also a steady state (P, N * ) = (0, 0) which represents the inflammation steady state (septic death state); this state will be analyzed in the next subsection by performing numerical simulations.

Linearization of system (4.9)-(4.10) in the neighborhood of the trivial steady state produces the system

. P = k pg P - k pm s m µ m P τ 1 , . N * = -µ n N * + s nr k np µ nr P τ 2 + s nr k nn µ nr N * τ 2 .
So the associated characteristic equation is given by

k pg -λ - k pm s m µ m e -λτ 1 0 s nr k np µ nr e -λτ 2 -µ n -λ + s nr k nn µ nr e -λτ 2 = 0, namely D(λ, τ 1 , τ 2 ) ≡ D 1 (λ, τ 1 ) • D 2 (λ, τ 2 ) = 0, (4.11)
where

D 1 (λ, τ 1 ) = λ -k pg + k pm s m µ m e -λτ 1 , D 2 (λ, τ 2 ) = λ + µ n - s nr k nn µ nr e -λτ 2 .
It is known from [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] that, when τ 1 = τ 2 = 0, the corresponding eigenvalues of (4.11) are real and given by λ = (k pg µ m -k pm s m ) /µ m and λ = (s nr k nn -µ n µ nr ) /µ nr . In what follows we assume that the inequality s nr k nn < µ n µ nr holds true. Consequently, in absence of delays, the steady state (0, 0) of system (4.9)-(4.10) is locally asymptotically stable if k pg < k pm s m /µ m .

Letting the time delays τ 1 and τ 2 varied, the trivial fixed point of system (4.9)-(4.10) may lose its stability. In order to consider the effects of the time delay, we need to investigate the boundary of the stability region determined by the equations λ = 0 and λ = iω (ω > 0). Letting λ = 0 in (4.11), one has that D(0, τ 1 , τ 2 ) = 0 having assumed that k pg µ m -k pm s m < 0. Thus, only the case λ = iω (ω > 0) needs to be analyzed.

Case τ 1 > 0 and τ 2 = 0 Eq. (4.11) becomes D(λ, τ 1 , 0) ≡ D 1 (λ, τ 1 ) • D 2 (λ, 0) = 0. Let λ = iω (ω > 0) be a root of D(λ, τ 1 , 0) = 0. Since D 2 (iω, 0) = 0, λ = iω has to solve D 1 (λ, τ 1 ) ≡ λ-k pg +(k pm s m /µ m ) e -λτ 1 = 0, i.e. D 1 (iω, τ 1 ) = 0. Setting a = k pg , b = -k pm s m /
µ m and τ = τ 1 , we note that this equation writes as λ -a -be -λτ = 0, which is Eq. (4.5). According to the analysis performed in the previous section, the following result holds. Proposition 4.3.1 Let τ 1 > 0, τ 2 = 0 and k pg < k pm s m /µ m . The steady state (P, N * ) = (0, 0) of system (4.9)-(4.10) is locally asymptotically stable for τ 1 < τ 1 0 and unstable for τ 1 > τ 1 0 , where

τ 1 0 = 1 ω 0 tan -1 ω 0 k pg and ω 0 = k pm s m µ m 2 -k 2 pg .
If τ 1 = τ 1 0 then the system (4.9)-(4.10) undergoes a Hopf bifurcation at (P, N * ) = (0, 0).

Case τ 1 = 0 and τ 2 > 0

In this case, Eq. (4.11) is D(λ, 0,

τ 2 ) ≡ D 1 (λ, 0) • D 2 (λ, τ 2 ) = 0. If λ = iω (ω > 0
) is a root of D(λ, 0, τ 2 ) = 0, then, being D 1 (iω, 0) = 0, we conclude that λ = iω is a solution of D 2 (λ, τ 2 ) ≡ λ + µ n -(s nr k nn /µ nr ) e -λτ 2 = 0, i.e. D 2 (iω, τ 2 ) = 0. Again, we can use Eq. (4.5) but with a = -µ n , b = s nr k nn /µ nr and τ = τ 2 . Contrary to the previous case, we find that the condition |b| > |a| is not satisfied. The fact that there is not purely imaginary root satisfying D(iω, 0, τ 2 ) = 0 leads to the following result.

Proposition 4.3.2 Let τ 1 = 0, τ 2 > 0. Then the steady state (P, N * ) = (0, 0) of system (4.9)-(4.10) is locally asymptotically stable.

Case τ 1 > 0 and τ 2 fixed in its stable interval

We now consider Eq. (4.11) with τ 2 in its interval of stability, regarding τ 1 as a parameter. Let λ = iω (ω > 0) be a root of (4.11). Then, D 1 (iω, τ 1 ) • D 2 (iω, τ 2 ) = 0. From the previous subsection, one has that D 2 (iω, τ 2 ) = 0, so that we must have D 1 (iω, τ 1 ) = 0. It is now straighforward that previous arguments imply the following result.

Proposition 4.3.3 Let τ 1 > 0, τ 2 > 0 and k pg < k pm s m /µ m . The steady state (P, N * ) = (0, 0) of system (4.9)-(4.10) is locally asymptotically stable for τ 1 < τ 1 0 and unstable for τ 1 > τ 1 0 , where

τ 1 0 = 1 ω 0 tan -1 ω 0 k pg and ω 0 = k pm s m µ m 2 -k 2 pg .
If τ 1 = τ 1 0 then system (4.9)-(4.10) undergoes a Hopf bifurcation at (P, N * ) = (0, 0).

Case τ 1 fixed in its stable interval and τ 2 > 0 This case states that τ 1 is in its stable interval [0, τ 1 0 ) and τ 2 is regarded as a parameter. Assume that Eq. (4.11) has purely imaginary solution of the form λ = iω (ω > 0). Then it can be seen that D 1 (iω, τ 1 ) = 0 since τ 1 ∈ [0, τ 1 0 ) and D 2 (iω, τ 2 ) = 0. Consequently, we have the following result.

Proposition 4.3.4 Let τ 1 ∈ [0, τ 1 0 ), τ 2 > 0 and k pg < k pm s m /µ m . Then the steady state (P, N * ) = (0, 0) of system (4.9)-(4.10) is locally asymptotically stable for τ 2 > 0.

Case

τ 1 = τ 2 = τ Eq. (4.11) becomes D(λ, τ ) ≡ D 1 (λ, τ ) • D 2 (λ, τ ) = 0. Once again, if λ = iω (ω > 0) is such that D(iω, τ ) = 0, we have that D 1 (iω, τ ) = 0.
Applying the analysis of previous subsections, we obtain conditions under which a family of periodic solutions bifurcate from the trivial equilibrium.

τ 0 = 1 ω 0 tan -1 ω 0 k pg and ω 0 = k pm s m µ m 2 -k 2 pg .
such that the steady state (P, N * ) = (0, 0) of system (4.9)-(4.10) is locally asymptotically stable for τ < τ 0 and unstable for τ > τ 0 . Moreover, for τ = τ 0 , the Hopf bifurcation occurs at (P, N * ) = (0, 0).

Numerical investigations

This section is devoted to further investigations on the model (4.9)-(4.10). Specifically by employing numerical solutions we perform the stability analysis and we investigate on the existence of Hopf bifurcation in the other fixed point (P, N * ) of the model (4.9)-(4.10) where analytical results have not been reported. Accordingly we perform a sensitivity analysis on the time delays τ 1 and τ 2 . The stability analysis is performed by considering a specific model and in order to compare our model with the model proposed in [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF], we set the parameters of the delayed model (4.9)-(4.10) as follows: The first step is to compute the fixed points. As known, the steady states are points where the nullclines intersect. As Figure 4.1 shows we have two different steady states: the fixed point E 0 = (0, 0) whose stability analysis has been performed in the previous section and the fixed point E 1 = (5 × 10 5 , 1.6). In what follows the numerical investigations focus on the stability analysis of the fixed point E 1 .

p ∞ =
The first set of simulations refers to the case τ 1 = τ 2 = τ . Letting τ varying, we have found that a Hopf bifurcation occurs at τ = 0.83358. Indeed as Figure 4.2 shows, if τ < τ the fixed point E 1 is asymptotically stable and if τ > τ the fixed point E 1 is unstable (the numerical solutions has been obtained for a set of initial conditions near the fixed point E 1 ). It is worth stressing that in the previous section we have proved analytically that the fixed point E 0 undergoes a Hopf bifurcation. This set of simulations is susceptible of biological interpretation. Indeed in the equations of P and N * at time t we have introduced also the role of these cells at time t -τ in order to take into account different stages of activation. Differently from [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF], our results show how it is important during the interactions to take care of the time at which the cells are activated. According to Figure 4.2, if at time t we have also cells whose activation stage is that at time t-τ , with τ < τ , the system will reach the septic death; if at time t we consider cells whose activation stage is that at time t -τ , with τ > τ , the system can reach the health stage (remember that in this case also the steady state E 0 undergoes a Hopf bifurcation).

The second set of numerical investigations refers to the case τ 1 = 0 and τ 2 ≥ 0. As Figure 4.3 shows, for a set of initial conditions near the fixed point E 1 , when τ 2 < τ the fixed point E 1 is always locally asymptotically stable; for τ 2 > τ the fixed point E 1 is unstable (in particular the origin is always asymptotically stable as already proved in the previous section). This set of results shows again how the time delay influences the asymptotic behavior of the system: health state or septic death. In the third case, namely when τ 1 ≥ 0 and τ 2 = 0, our numerical results show that the fixed point E 1 is always asymptotically stable. The last set of numerical simulations refers to the sensitivity analysis with respect the parameter k pg and in the case τ 1 = τ 2 = τ . In the case τ = 0.4, the Figure 4.4 shows that for k pg = 2.7 the system reaches the steady state E 0 and for k pg = 2.95 the system reaches the steady state E 1 .

Conclusions and future perspective

The present chapter has been devoted to generalize the mathematical model developed in [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation[END_REF] by inserting two time delays in order to take into account cells with different stage of activation.

The asymptotic analysis has been concerned with the stability analysis of the steady states and the sufficient conditions under which a Hopf bifurcation occurs. The analysis has shown how time delays can influence the whole dynamics and in particular the stability of the steady state.

In particular the numerical results and the bifurcation analysis suggest that the magnitude of the time delay plays several important roles in the restoration of health.

It is worth stressing that the biological relevance of the analysis performed in the present chapter and the related conclusions are limited by the simplifications present in the proposed delayed model. Therefore the development of a much larger model that is also able to reproduce experimental data is object of future investigations.

Research directions include the derivation of an explicit algorithm for determining the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions. This step can be pursued by employing the center manifold theory and the normal form method [START_REF] Hassard | Theory and Applications of Hopf Bifurcation[END_REF]. Further research directions include generalizations of the delayed model proposed in the present chapter, which takes into account tissue damage and anti-inflammatory mediators such as cortisol and interleukin-10. Moreover the introduction of the adaptive branch of the immune system is object of further investigations.

Finally the derivation of a mathematical model for the acute inflammatory response to infection based on the thermostated kinetic theory [START_REF] Bianca | Existence of stationary solutions in kinetic models with gaussian thermostats[END_REF][START_REF] Bianca | High-order moments conservation in thermostatted kinetic models[END_REF] is part of future research perspective.

Mimic Therapeutic Actions against Keloid by thermostated Kinetic Theory Methods

This chapter deals with the modeling of a wound healing disease under a therapeutic action by employing the methods of the thermostated kinetic theory for active particles. In particular, in order to test a therapeutic action for the keloid formation and the possible development of a cancer, an external force field coupled to a Gaussian thermostat is introduced into a mathematical model recently proposed. Specifically the model depicts the competition of the immune system cells with a virus, the mutated fibroblast cells, and the cancer cells. Employing a computational analysis, the effect of three different external force fields mimic therapeutic actions is analyzed: A vaccine for the virus, the PUVA therapy for the keloid and a vaccine for the cancer. The results are in agreement with the evidence that the sole action of the immune system is not sufficient to obtain a total depletion of keloid. Further refinements and developments are also discussed into the chapter. It is worth pointing out that the contents of this chapter constitutes the paper [START_REF] Bianca | Mimic therapeutic actions against keloid by thermostatted kinetic theory methods, to appear[END_REF].

Introduction

As already mentioned in chapter 5, wound healing is a dynamic process consisting of four continuous, overlapping, and precisely programmed phases: Hemostasis, the inflammation phase, the proliferation phase, and the remodelling phase, see the review paper [START_REF] Guo | Factors affecting wound healing[END_REF]. The events of each phase must happen in a precise and regulated manner. Interruptions, aberrancies, or prolongation in the process can lead to delayed wound healing or a non-healing chronic wound, such as keloid. Keloid is a hyperproliferative response of connective tissue in response to skin trauma. The causes which trigger this phenomenon are poorly understood and currently no successful treatment have been developed [START_REF] Kelly | Medical and surgical therapies for keloids[END_REF]. In particular the defective control mechanisms in keloid result in the excessive cell proliferation and extracellular matrix synthesis: keloidderived fibroblasts have a greater proliferative capacity than normal derived fibroblasts [START_REF] Calderon | Increased proliferation in keloid fibroblasts wounded in vitro[END_REF]. Even if fibroblasts have a major role, other cells like keratinocytes and melanocytes can be involved [START_REF] Kelly | Medical and surgical therapies for keloids[END_REF], and the causes can be also the presence of a virus [START_REF] Alonso | Keloids: A viral hypothesis[END_REF] and a generic susceptibility [START_REF] Saed | Analysis of p53 gene mutations in keloids using polymerase chain reaction-based single-strand conformational polymorphism and DNA sequencing[END_REF]. The possible therapy for keloid strictly depends on the location, size, depth of the lesion and the age of the patient. Therapeutic treatment includes occlusive dressings, compression therapy, intralesional corticosteroid injections, cryosurgery, excision, radiation therapy, laser therapy, interferon therapy, 5-fluorouracil, doxorubicin, bleomycin, verapamil, retinoic acid. Other promising therapies include antiangiogenic factors, including vascular endothelial growth factor inhibitors, phototherapy, tumor necrosis factor (TNF)-alpha inhibitors, and recombinant human interleukin, which are directed at decreasing collagen synthesis.

Recently the defective induction of stress-induced premature senescence during wound healing has been proposed as a possible mechanism of keloid formation [START_REF] Blazic | Defective induction of senescence during wound healing is a possible mechanism of keloid formation[END_REF]. Specifically keloid fibroblasts undergo senescence at a rate lower than that of normal scar, thus depositing collagen and other extracellular matrix proteins beyond that expected in normal wound healing. The proposed mechanism could lead to new treatment possibilities for keloid, e.g. a therapy that induce senescence could be used to prevent the formation of keloid, and consecutively enable the formation of a normal scar.

From the modeling point of view, the fibro-proliferative disease has been widely studied at macroscopic level, but insight at cellular and sub-cellular level scale is still lacking, see the review paper [START_REF] Fusi | Macroscopi models for fibroproliferative disorders: A review[END_REF] and references cited therein. Therefore, there is a remarkable need of developing the necessary modeling approach that can handle the full range of time and length scales required to model complex biological systems, with emphasis on spatially distributed systems, and apply this framework to fibrosis disease such as keloid.

The keloid formation has been recently modeled in [START_REF] Bianca | Mathematical modelling for keloid formation triggered by virus: Malignant effects and immune system competition[END_REF] and further analyzed in [START_REF] Bianca | Bifurcation diagrams for the moments of a kinetic type model of keloid -immune system competition[END_REF]. The model, which refers to the keloid formation triggered by a virus and its possible malignant effects under the immune system surveillance, is based on the kinetic theory for active particles [START_REF] Bianca | Towards a Mathematical Theory of Multiscale Complex Biological Systems[END_REF]. At the basis of the methods there is the decomposition of the whole system into different functional subsystems composed by particles, called active particles, able to express a function or strategy. Therefore the microscopic state of the active particles includes, in addition to the classical geometrical and mechanical variables, also the activity variable. The description of each functional subsystem is obtained by the definition of a distribution function over the microscopic state of interacting particles. The evolution of the system is determined by microscopic nonlinear interactions, which take into account not only mechanical rules but also modifications of the biological state. It is worth stressing that the kinetic theory for active particles has been also employed for the modeling of other biological systems, see, among others, papers [START_REF] Jorcyk | Experimental versus numerical data for breast cancer progression[END_REF][START_REF] Kolev | Adaptive immunity and CLT differentiationa kinetic modeling approach[END_REF][START_REF] Kolev | Numerical simulations for tumor and cellular immune system interactions in lung cancer treatment[END_REF][START_REF] Cattani | A Class of Solutions for the Hybrid Kinetic Model in the Tumor-Immune System Competition[END_REF].

This chapter is concerned with a further refinement of the mathematical model proposed in [START_REF] Bianca | Mathematical modelling for keloid formation triggered by virus: Malignant effects and immune system competition[END_REF]. Specifically in order to take into account possible therapeutical actions, the model [START_REF] Bianca | Mathematical modelling for keloid formation triggered by virus: Malignant effects and immune system competition[END_REF] is generalized by introducing different external force fields, which mimic a therapy at the macroscopic scale. The introduction of an external force field moves the system out of equilibrium and in order to have a steady state a Gaussian thermostat is coupled to the external force field. The Gaussian thermostat is thus a damping term adjusted to control the energy into the system, the interested reader in a more deeper understanding of the thermostat and its applications in molecular dynamics simulations is referred to the book [START_REF] Evans | Statistical mechanics of nonequilibrium fluids[END_REF] and to the recent review paper [START_REF] Jepps | Deterministic thermostats[END_REF]. The resulting framework is thus a thermostated kinetic model where the excess of energy introduced by the external force field is removed by the thermostat. This new framework acts as a general paradigm for the derivation of a specific model for keloid. In particular the thermostated kinetic model proposed in the present chapter is the object of a computational analysis to test a recent hypothesis of therapeutic action for the keloid formation and the possible development of a cancer. The model depicts the competition of the immune system cells with a virus, the mutated fibroblast cells, and the cancer cells by means of the density function and the distribution function of cells. Three different external forces, mimic therapeutic actions, are taken into account: a vaccine for the virus, the PUVA therapy for the keloid and a vaccine for the cancer cells. The results are in agreement with the evidence that the sole action of the immune system is not sufficient to obtain a total depletion of keloid.

The contents of the present chapter are outlined as follows: After this introduction, Section 7.2 is concerned with the presentation of the thermostated kinetic theory for active particles and specifically with the derivation of a mathematical structure that acts as a general paradigm for the derivation of specific models. The framework includes an external force field coupled with a Gaussian thermostat that allows the control of the global activation energy. Section 7.3, after a short phenomenological analysis with the aim of understanding who are the main actors that are responsible for the keloid formation, is devoted to the derivation of our model for keloid formation and specifically to the definition of the functional subsystems, the interactions among the cells and the introduction of the external force fields. Section 7.4 deals with a computational analysis on a test case obtained by setting a value to the different parameters of the model. The computational analysis consists of a sensitivity analysis on the magnitude of the external force fields when the rate of heterogeneity of cells is let vary. Finally Section 7.5 concludes the chapter by suggesting possible further refinements and developments, to be regarded as future research directions.

The underlying thermostated kinetic framework

This section is devoted the derivation of the thermostated kinetic framework that can be proposed for the modelling of most real-world complex systems. The section is presented through two sequential subsections. First we deal with the complexity problem, the decomposition of system into subsystems and their statistical representation, and finally with the derivation of the relative class of evolution equations that act as a general paradigm for the derivation of specific models.

Complexity, functional subsystems decomposition and representation

Complexity is an intrinsic characteristic of most living systems that makes the modelling of the system disputable. Indeed not only the large number of elements constituting the system but also the emergence of behaviour, that arises as result of the whole interactions among the elements that occur in nonlinear fashion, are key issues of the complexity. The complex behavior is also due to the fact that the living entities (cells, animals, pedestrians), differently from the inert matter entities (billard balls, gas particles, electrons), have the ability to perform specific strategies and functions such that small variations in their will, could modify the overall asymptotic dynamics. Moreover external factors influence their strategies and consequently the occurrence of their emerging behaviour. Finally the strategy of entities at the lower scale (microscopic or mesoscopic) determines the behaviors of the system at the upper scale (macroscopic). From the modelling viewpoint, the onset of complexity indicates the needing to develop tools and methods which allow to reduce the complexity.

In the biological systems case, several authors suggest the approach of system biology. In what follows we use the systems approach proposed in [START_REF] Bianca | Towards a Mathematical Theory of Multiscale Complex Biological Systems[END_REF]. Accordingly the overall system S is viewed as composed by a large number of elements (particles of the inert matter, animals, cells, pedestrians) that interact in nonlinear matter and under the effect of macroscopic external force fields. In this approach the elements of the system are called active particles because in addition to their possibility to be allocated, at time t, in the position x with velocity v, they also have an active role in the evolution of the system (the ability to express strategy) modelled by the scalar variable u ∈ D u . The triplet (x, v, u) is called the microscopic state of the particles. The active particles are of heterogeneous type and they are selected for the ability they express and not for their size and shape. Particles having the same ability are grouped into a subsystem, called functional subsystem. Therefore the overall complex system S is decomposed into different functional subsystems S i , for i ∈ {1, 2, . . . , n}, such that: S i ∩ S j = ∅, ∀i = j, and S = n i=1 S i .

(5.1)

It is worth precising that the decomposition into functional subsystems is a flexible approach to be adapted at each system and is related to the activity variable. Finally a functional subsystem could be itself a complex system and therefore can be decomposed into subsystems. The evolution of each functional subsystem S i is represented by the statistical mechanics approach. Accordingly, we assume that the active particle during a time interval [0, T ], or [0, +∞), is located in x which attains values in a (bounded or unbounded) domain D x ⊂ R d (usually d = 3), with velocity variable v that attains values in a domain D v ⊂ R d and activity variable u which varies in a domain D u ⊂ R. Then the kinetic of the ith functional subsystem, for i ∈ {1, 2, . . . , n}, is depicted by the evolution of the distribution function

f i (t, x, v, u) defined on [0, T ] × D x × D v × D u .
For any fixed time t, the quantity f i (t, x, v, u)dx dv du stands for the density of particles in the volume element dx dv du centered at (x, v, u). Let Ω = D x × D v × D u be the domain of all the possible microscopic states and dΩ = dx dv du the Lebesgue measure on Ω, then the minimal assumption that one can make on f is that for all t 0,

f i (t, •, •, •) ∈ L 1 loc (Ω, dΩ).

The homogeneous thermostated kinetic framework

This section is concerned with the underlying thermostated kinetic framework which acts as a general paradigm for the derivation of specific models. Specifically a large number of homogeneously distributed cells (active particles) interact in nonlinear matter and macroscopic external force fields F i = F i (u) : D u → R + extent an action on the system. The microscopic state of a cell is the variable u, which means that the cell at time t ∈ [0, ∞), has activity variable u ∈ D u ⊂ R. Cells expressing the same function are grouped into a subsystem, called functional subsystem. The evolution of each functional subsystem is depicted by the distribution function f i (t, u) : [0, ∞)×D u → R + , for i ∈ {1, 2, . . . , n}, and such that, for any fixed time t, the quantity f i (t, u)du represents the density of cells in the elementary volume du centered at u. Let f = f (t, u) = (f 1 (t, u), . . . , f n (t, u)) be the vector whose components are the distribution functions of the functional subsystems and

f (t, u) = n i=1 f i (t, u).
(5.

2)

The evolution equation for the ith functional subsystem is obtained by equating the time derivative of f i to the balance of the inlet and outlet flows in the elementary volume du. Accordingly we have:

∂ t f i + ∂ u F i 1 -u Du u f du f i = J i [f ], (5.3)
where

J i [f ] = J i [f ](t, u)
is the operator that models the gain-loss of cells due to transitions in the activity variable (conservative term) and it reads:

J i [f ] = n j=1 Du×Du η ij A ij (u * , u * , u)f i (t, u * ) f j (t, u * ) du * du * -f i (t, u) n j=1 Du η ij f j (t, u * ) du * .
(5.4)

where:

• η ij models the probability that a cell of the ith functional subsystem with activity u * interacts instantaneously with a cell of the jth functional subsystem with activity u * ;

• A ij = A ij (u * , u * , u) : D u × D u × D u → R +
is the density function modeling the probability that cells of the ith functional subsystem, with activity u * , interacting with cells of the jth functional subsystem, with activity u * , reach the activity u. In particular A ij (u * , u * , u) satisfies the following identity: ; the possibility of this interaction can be measure by introducing the non-negative function η ij (u * , u * ) which represents the interaction rate; the probability that after this interaction the candidate cell undergoes a change in its microscopic state is measured by the non-negative function A ij (u * , u * , u). Bearing all above in mind, the (infinitesimal) result of the interaction reads:

Du A ij (u * , u * , u) du = 1, ∀ u * , u * ∈ D u . Remark
A ij (u * , u * , u) η ij (u * , u * ) f i (t, u * ) du * f j (t, u * ) du * , (5.6)
and summing up with respect to all the cells with activities u * and u * , we obtain the operator which models the gain of cells with activity u, see Eq (5.4).

The operator T

F i [f ] = T F i [f ](t, u)
is the transport term that models the Gaussian thermostat [129,[START_REF] Ruelle | Smooth Dynamics and New Theoretical Ideas in Nonequilibrium Statistical Mechanics[END_REF], and it reads:

T F i [f ] := ∂ u F i 1 -u Du u f (t, u) du f i (t, u) . (5.7)
In particular (5.7) is a damping operator adjusted to control the global activation energy:

E 2 [f ](t) = n i=1 Du u 2 f i (t, u) du.
The framework (5.3) can be further generalized by introducing the role of nonconservative processes. Specifically, interactions among the cells may generate proliferation/destruction of other cells (birth-death process). This type of interaction is modeled by the following operator

N i [f ] = N [f ](t, u): N i [f ] = f i (t, u) n j=1 Du α ij f j (t, u * ) du * , (5.8)
where α ij = η ij µ ij , being µ ij the net proliferative/destructive rate.

Because of DNA corruptions, cells can become cells of another functional subsystem. These kinds of interactions are modeled by the following operator

M i [f ] = M i [f ](t, u): M i [f ] = n h=1 n k=1 Du×Du β ihk f h (t, u * ) f k (t, u * ) du * du * .
(5.9)

where β ihk = η hk ϕ i hk , being ϕ i hk the net mutative rate into the ith functional subsystem, due to interactions that occur with rate η hk between the cells with activity u * of the hth functional subsystem and the cells with activity u * of the kth functional subsystem. Bearing all above in mind, the thermostated kinetic framework with proliferative/destructive and mutative interactions reads:

∂ t f i (t, u) + T F i [f ](t, u) = J i [f ](t, u) + N i [f ](t, u) + M i [f ](t, u).
(5.10)

It is worth stressing that the parameters of the model can be function of the activity, namely

α ij (u * , u * ) = η ij (u * , u * ) µ ij (u * , u * ) and β ihk (u * , u * ) = η hk (u * , u * ) ϕ i hk (u * , u * ).
However in order to simplify our model the parameters α ij and β ihk will be assumed as constants.

From the mathematical point of view, the Cauchy problem related to the general framework (5.10) has been analyzed and the existence and uniqueness of the solution has been proved including the existence of stationary solutions [START_REF] Bianca | High-order moments conservation in thermostatted kinetic models[END_REF]. Therefore we are allowed to develop the appropriate computational methods to obtain simulations of a specific model.

A thermostated kinetic model for keloid with therapy

This section, according to the phenomenological description of [START_REF] Bianca | Mathematical modelling for keloid formation triggered by virus: Malignant effects and immune system competition[END_REF], deals with the derivation, at the cellular scale, of a model for the keloid formation with therapy in the framework of the thermostated kinetic theory. According to the general framework (5.10), first we deal with the characterization of the functional subsystems and their statistical representation, then with the modelling of the microscopic interactions among the selected functional subsystems, and finally with the derivation of the relative class of evolution equations.

Functional Subsystems and their Activity Variable

The model is developed assuming that the keloid formation involves four interacting functional subsystems and these interactions may also generate a subsystem constituted by cancer cells. Each subsystem is able to express a well-defined biological function represented by the scalar variable u ∈ D u = [0, ∞[. Specifically we distinguish among the following subsystems:

1 -Fibroblast cells (Fc). The activity variable u represents the proliferation ability. A fibroblast cell may undergo a mutation as consequence of a viral action. This mutation gives the cells a significant advantage with regard to proliferation, and allows it and its descendant to quickly advance along the keloid. Thus normal fibroblasts cells may generate cells in a new subsystem that we call, according to medical terminology, keloid fibroblast cells.

-Viruses (V).

The activity variable u is a magnitude of their aggressiveness related to their proliferation ability. Virus, as known, is an infectious agent that can replicate only inside the cells which, as consequence of the infections, became viral cells or activated viruses. Virus infections eventually result in the death of the host cell. When the adaptive immune system of a vertebrate encounters a virus, it produces specific antibodies that bind to the virus and render it non-infectious (humoral immunity). The host cell is destroyed by killer T cells (cells-mediated immunity).

-Keloid-Fibroblast cells (KFc).

The activity variable u represents the proliferation ability but these cells are characterized by a relatively greater ability to proliferate with respect to the normal fibroblast cells (genetic instability). Increasing values of u indicate an increasing intensity of proliferation.

-Cancer cells (Cc).

The activity variable u is a magnitude of their progression ability. These cells are characterized by a greater and greater genetic instability with respect to the KFc. Increasing values of u take into account these capabilities.

-Immune system cells (ISc).

The activity variable u is a magnitude of the activation and of the response for the immune system to foreign agents. These cells have the ability to contrast KFc, Cc, and V through recognition and destruction.

In the sequel, when we refer to the subsystems 2, 3, and 4 as a whole, we will call them non-self cells. In table 5.1 we have summarized the funtional subsystems, the related activity variables, and the distribution functions of the cells involved in the model.

Functional Subsystems Activity Distribution Function

Fibroblasts cells (Fc) Proliferation f 1 (t, u)

Viruses (V) Aggressiveness f 2 (t, u) Keloid Fibroblasts cells (KFc) Proliferation f 3 (t, u) Cancer cells (Cc) Progression f 4 (t, u)
Immune System cells (ISc) Activation f 5 (t, u)

Table 5.1: The functional subsystems, the activity variables, and distribution functions of the cells involved in the model.

The Microscopic Interaction Terms

The derivation of evolution equations for the distribution function of each subsystem needs the modelling of the microscopic interactions among individuals of the various subsystems. Interactions which have non trivial outputs, i.e. an effective change (either in the microscopic state or in the entities number) occurs after the interaction, are not considered. Thus biological phenomena corresponding to interactions that give trivial outputs will not be taken into account, as for instance the physiological birth and death of the entities. The microscopic interaction terms are derived under the following mathematical assumptions:

1 -The overall state of the functional subsystems is described by the distribution functions f 1 (t, u) (fibroblast cells), f 2 (t, u) (viruses), f 3 (t, u) (keloid fibroblast cells), f 4 (t, u) (cancer cells), and f 5 (t, u) (immune system cells).

-

The local density and the local activation for each subsystem i ∈ {1, 2, 3, 4, 5} write respectively:

E 0 [f i ](t) = ∞ 0 f i (t, u) du, E 1 [f i ](t) = ∞ 0 uf i (t, u) du . (5.11)
destructive interactions, the parameters γ and λ are related to transition encounters, and ε is a scale parameter, see Table 5.2. The all parameters have to be regarded as positive constants (eventually equal to zero), small with respect to unity, and have to be tuned by suitable experiments.

Interactions Fc V KFc Cc ISc Fc Proliferative ε 2 β Destructive -εδ Mutative εγ γ V Conservative εα Proliferative εβ εβ Destructive -δ KFc Conservative α Proliferative β β Destructive -ε 2 δ -ε 2 δ Mutative λ Cc Conservative ε 2 α Proliferative εβ Destructive -δ ISc Proliferative β i ε 2 β i β i Destructive δ i -ε 2 δ i δ i
Table 5.2: Functional subsystems, interaction terms, and related parameters.

Computational analysis: Mimic therapeutic actions

This section is concerned with the computational analysis of the model (5.21) and specifically with the evolution of the different functional subsystems when an external action acts on the system as a therapeutic action. The main aim is to simulate the prompt response against the formation and evolution of keloid and the possible onset of cancer. The computational analysis is thus addressed to analyze the effects of three different therapeutic actions: An action which mimics a vaccine against the virus, an action which mimics a vaccine against cancer cells, and an action which can also mimic surgery on keloid. It is important to note that, according to our model, the onset of cancer cells is a consequence of mutations in the keloid cells because of the virus, the latter in part also responsible for mutations in the fibroblasts cells which generate the keloid (remember that according to our assumptions, the genetic susceptibility is also responsible for keloid formation). Therefore the main role of the external actions should be to act on the functional subsystems of virus and keloid cells. However in order to have a more global view of what the model developed in the present chapter is able to reproduce, we will consider also an external action on the cancer cells.

It is worth stressing that our model is an exploratory model. Thus, at this stage, we are only interested in the emerging phenomena that the model is able to reproduce and not in the tuning with an in-vivo or in-vitro experiment or empirical data, which can be considered as subject of further investigations. Accordingly, the computational analysis focuses on the model (5.21) when only one external action is applied. The dynamics of the model when the three actions can act at the same moment will be straightforward. However to think that the three different actions may be applied at the same time can not be suitable then we believe that our model will fit well if a combination of the three different actions will be performed at different times. Specifically the first step is to act against the virus (F 2 = 0, F 3 = F 4 = 0); when the virus is killed, the reached nonequilibrium stationary state will be the initial state for the model (5.21) with the action of a therapy for keloid (F 3 = 0, F 4 = 0); when keloid cells are eliminated, the new reached nonequilibrium stationary state will be the initial state for the model (5.21) with the therapeutic action against the cancer cells (F 4 = 0). In particular the intermediate step can be performed thanks to the introduction of the thermostat that allows the existence of a nonequilibrium stationary state.

The computational analysis will be addressed by fixing nine of the eleven parameters and performing a sensitivity analysis on the parameter α and on the therapeutical action F i , for i ∈ {2, 3, 4}. In the analysis we will show the evolution of the density, the behavior of the distribution function and the activation energy. The computational scheme is that of the wellknown generalized collocation method where the variable u is discretized into a suitable set of collocation points. The integral terms are approximated by means of algebraic weighted sums in the nodal points of the discretization. The particularization of the evolution equations in each node and the enforcing of the initial conditions transform the model that is a system of integro-differential equations into a systems of ordinary differential equations, describing the evolution of the values of the distribution functions in the node of the collocation, see Section 2 of chapter [START_REF] Bianca | Bifurcation diagrams for the moments of a kinetic type model of keloid -immune system competition[END_REF] for the all details.

The choice of the distribution functions at time t = 0 is based on the assumption that, before the formation of keloid and the related onset of cancer, the virus infects the Fc. Our analysis starts when the number of Fc in the wound is equal to the number of V and a number of immune system cells have reached the wound. Accordingly we assume nonzero initial conditions only for the functional subsystem of Fc, V, and ISc.

The test case is based on the choice of the parameters as follows: γ = 0.4 (the rate of mutation of Fc into KFc is not negligible), δ = 0.3 (the ability of the immune system cells to inhibit the non-self cells is quite low), δ i = 0.5 (the non-self cells have an intermediate ability to inhibit the response of the immune system), β = 0.4 (the non-self cells have an intermediate ability to proliferate), β i = 0.35 (the rate of proliferation of ISc is quite low), λ = 0.5 (the rate of mutation of KFc into Cc is not negligible) and = 0.5.

It is worth stressing that the computational analysis performed in the present chapter do not cover the whole variety of conceivable dynamics but represents a useful test case.

Simulating the effects of a vaccine for the virus

This subsection deals with the computational analysis for the model (5.21) when a constant external force, mimic a vaccine for the virus, is introduced. Following the interaction rules proposed into the model, the onset of keloid is a consequence of the virus, then it is expected that if the vaccine is able to reduce the action of the virus then keloid formation and cancer can be prevented. Accordingly we have F 2 = 0, F 3 = F 4 = 0 and we let the magnitude of the parameter α vary from low to higher values, namely from low to higher heterogeneity. The results of the computational analysis are summarized as follows.

Analysis for low values of α. The effect of the external action is evident by comparing the dynamics depicted by the model (5.21) when F 2 = 0 with the dynamics when F 2 = 0.

The dynamics of the model (5.21) for F 2 = F 3 = F 4 = 0 has been widely debated in [START_REF] Bianca | Mathematical modelling for keloid formation triggered by virus: Malignant effects and immune system competition[END_REF], therefore the new computational analysis will be focused on the case F 2 = 0 and F 3 = F 4 = 0, and more precisely F 2 ∈ {0, 0.0005, 0.01}. Looking at the left panel of Figure 5.1, for α = 0.2 (low rate of heterogeneity) the heterogeneity of virus is bounded, then the main effect of the external force (vaccine for the virus) is to decrease the maximum of the density of V when the magnitude of the vaccine increases. Therefore, thanks also to the action of the immune system, we can see a faster depletion of V and consequently of KFc and Cc. However, as the right panel of Figure 5.1 shows, the keloid cells start to increase again; this is a consequence of the genetic susceptibility that thus requires the definition of a therapeutical action for keloid. The global energy of the system even in the presence of the external action is preserved thanks to the thermostat. Analysis for intermediate values of α. Setting α = 0.5, the heterogeneity of the non-self cells becomes less negligible with respect the previous case. Therefore the role of the external force is now fundamental for avoiding keloid formation and onset of cancer. The computational analysis shows that high values of the vaccine are required in order to reduce the number density of keloid cells with high values of activity. Consequently the number density of Cc is reduced. The vaccine for the virus has an important role in the proliferation of the immune system cells. Indeed as Figure 5.2 shows, in absence of a vaccine for the V the immune system is inhibited; the vaccine helps the immune system cells to proliferate again. Analysis for high values of α. The heterogeneity of the non-self cells in now very high thus keloid formation and cancer need to be inspected. As Figure 5.3 shows the vaccine is able to control the evolution of the virus and the depletion of virus with high levels of aggressiveness depends on the magnitude of the vaccine. However high values of the vaccine are necessary to eliminate keloid and cancer, see Figure 5.4, then the definition of a therapeutic action for keloid and cancer is now necessary.

Simulating the effects of a therapy for the keloid

This subsection deals with the computational analysis for the model (5.21) when a constant external force, mimic a therapeutical action for the keloid, is considered. Following [START_REF] Blazic | Defective induction of senescence during wound healing is a possible mechanism of keloid formation[END_REF], a therapy that induces senescence could be used to prevent the formation of a normal scar. Specifically the treatment proposed in [START_REF] Blazic | Defective induction of senescence during wound healing is a possible mechanism of keloid formation[END_REF] is based on photodynamic and PUVA therapy, which are capable of inducing cell senescence. The therapy, based on a combination of a photosensitizer and a light source, produces oxidative stress and thus produce higher quantity of senescent cells with the minor apoptotic and necrotic effect.

Bearing all above in mind, we assume that F 3 = 0, F 2 = F 4 = 0 and again we let the magnitude of the parameter α vary from low to higher values. Specifically for low values of α, the therapy and the immune system are capable of performing a prompt action, see the left panel of Figure 5.5 where the density of the KFc is depicted for different values of the external force. Moreover, as the right panel of Figure 5.5 shows, the global activation energy is controlled by thermostated. The effects of the therapy are visualized in Figure 5.6, where the distribution function of the KFc is depicted in absence of therapy and in the presence of the therapy. As Figure 5.6 shows, the therapy allows the depletion of cells at different stages of mutation and a bounded magnitude of therapy is required in order to have the total depletion of keloid. Moreover, differently from the case of a vaccine for the virus, in this case the increase of new keloid cells, after the elimination of keloid, is not observed. For intermediate values of α, the general effect of the therapy is that of acting homogeneously on the KFc and Cc and thus to obtain a distribution more uniform with respect to the degree of mutation. In particular the Fc are able to proliferate again. For high values of α, the rate of mutation allows the onset of keloid cells with high proliferation rate and of cancer cells with high values of progression. In particular even if the therapy allows the immune system to have a prompt response (see Figure 5.7) against the cancer cells (and eventual virus), the immune system cells are not able to prevent tumor formation at tissue scale. Therefore a high magnitude of the therapy is required. This is the case where surgery is the only therapy that can work.

Simulating the effects of a vaccine for the cancer

This subsection deals with the computational analysis for the model (5.21) when a constant external force, mimic a therapeutical action for the cancer, acts on the system. The vaccine for the cancer strictly depends on the type of tumor developed. Therefore the main aim of this subsection is to show that our model is able to consider the introduction of a vaccine and to simulate the evolution of the cancer. In particular we will focus on the case of high values of α. According to our model, the development of cancer is the final result of the keloid formation, then the vaccine for cancer has a fundamental action only on the Cc and the ISc. As Figure 5.8 shows, the vaccine is able to inhibit the formation of a tumor at the macroscopic scale. Moreover the immune system is able to proliferate again, see the right panel of Figure 5.9. Finally, as the right panel of Figure 5.9 shows, the action of the therapy allows to maintain bounded the global activation energy of system (thanks to the introduction of the thermostat).

Conclusions and research perspectives

The present chapter has been devoted to test the capability of a new thermostated kinetic framework for the active particles to model a therapy for the keloid formation and the possible development of cancer. By employing different therapies on the different functional subsystems that compose the system, we have shown that some therapies proposed in the pertinent literature can be well described by our model. The analysis performed in this chapter has been of computational kind and has been focused on the reaching of the emerging phenomena that are typical of keloid formation. The role of the immune system has been taken into account as a whole system without specifying what are the cells involved in the competition; this assumption is based on the fact that its action on keloid is limited by the genetic susceptibility that does not allow to the immune system to recognize these cells as foreign cells and thus acts for a possible total depletion.

The model developed and analyzed in the present chapter is based on the introduction of an external force field, which mimics a therapy, at the macroscopic scale but in the absence of external interactions at the microscopic scale. The therapy acts directly on the number of foreign cells without performing an interaction at the microscopic level of the activity variable. The modeling of external therapy at the microscopic (cellular) scale can be performed by representing the external therapy as a functional subsystem that have the ability to modify the state u of the cell by a particular action related to the variable ω ∈ D u . Assuming that the ith inner functional subsystem interacts with the rth external agent, for r ∈ {1, 2, . . . , m}, and denoting by g ir = g ir (t, ω) : [0, ∞) × D u → R + the related distribution function (known function of its arguments), the microscopic external actions are modeled by the following operator where g i = (g i1 , . . . , g ir ) and

Q i [f , g i ] = Q i [f , g i ](t, u): Q i [f , g i ] = m r=1
• η e ir is the inner-outer encounter rate between the rth external agent, with state ω * , and the cell of the ith population, with state u * .

• B ir (u * , ω * , u) is the inner-outer transition probability density which describes the probability density that a cell of the ith population, with state u * , falls into the state u after an interaction with the rth external agent whose state is ω * . The density B ir satisfies, for all r ∈ {1, 2, . . . , m} and i ∈ {1, 2, . . . , n}, the following condition:

Du B ir (u * , ω * , u) du = 1, ∀ u * , ω * ∈ D u .

(5.23)

Bearing all above in mind, the thermostated kinetic framework for open systems with proliferative/destructive and mutative interactions now reads:

∂ t f i (t, u) + T F i [f ](t, u) = J i [f ](t, u) + N i [f ](t, u) + M i [f ](t, u) + Q i [f , g i ](t, u). (5.24)
It is worth stressing that usually the interaction domain of the cell with state u * is not the whole domain D u but a subset Ω u * ⊆ D u , which contains the cells with activity u * ∈ Ω u * that are able to interact with the particles with activity u * . This is a phenomenon that is typical in tumor dynamics and specifically when the immune system is not capable of interacting with the all tumor cells (tumor escape), see the review paper [START_REF] Malmberg | Escape from immune-and nonimmunemediated tumor surveillance[END_REF]. Specifically the immune response fails to completely eliminate the tumor, and the interact process results in the selection of tumor cell variants that are able to resist, avoid, or suppress the antitumor immune response, leading to the escape phase. Accordingly, a positive function ω(u * , u * ) can be introduced to weight the interactions among the cells; this function is assumed normalized with respect to integration over u * and it has a compact support in the domain of influence Ω u * ⊆ D u of the interactions. Moreover: (5.25)

The development of a model for the therapy for keloid that takes into account the above raised issues is object of future research directions. Moreover the possibility to perform an asymptotic analysis that allows the derivation of the dynamics at the tissue scale is a further research perspective. The asymptotic analysis can be performed by employing the methods developed in papers [START_REF] Degond | Mass and energy balance laws derived from high-field limits of thermostatted boltzmann equations[END_REF][START_REF] Bianca | Kinetic models coupled with gaussian thermostats: Macroscopic frameworks[END_REF]. Moreover the final goal is the identification of the parameters with the aim to compare the emerging phenomena with the experimental data.

Conclusion

This PhD thesis is concerned with general methods and models which can be applied to distinct scientific fields, as biology or finance. On the one hand, the scientist observes the data in distinct fields, and sees that there are analogous underlying characteritics, resulting in similar mathematical properties.

In the first part, the observed common phenomenon is algebraic: the degree vector and the principal eigenvector of a matrix are close componentwise, in both biological and financial data. From this observation, we analytically obtained an upper bound from the error vector, whose coordinates are the substraction (or division) of the principal eigenvector and the degree vector componentwise. This result has also been further tested for some known networks and financial data.

Using some of the mathematical results presented in the first part, we have extracted from a chromosomal contact map (i) the functional compartments in the human embryonic stem cell nucleus, and (ii) the 3D-structure of DNA in chromosomes.

In the second part, it has been shown that the acute inflammatory response is related to a bifurcation when delays are considered. We also mathematically described the mimic of the effect of a vaccine on a population of viruses.

Taken together, the work described in this thesis holds towards interesting perspectives. On the financial viewpoint, these perspectives involve trading strategies that will take advantage on the mathematical properties described in part one. Some further research perspectives would be to consider the random matrix characteristics linked to the closeness of the degree vector and the principal eigenvector. Some other research perspective deal with the comparison of the kinetic theory results with effective data.

Definition 1 . 1 . 2 11 )Remark 1 . 1 . 3

 11211113 Let M ∈ SIR n (R + ). We define Remark 1.1.2 Since M is irreducible, the numbers α and β are in R * + . Moreover n β (M e) (i) αn, ∀i ∈ {1, 2, . . . , n}. If αβ = 1 then max i∈{1,2,...,n} (M e) (i) = min i∈{1,2,...,n}

Figure 1 . 1 :

 11 Figure 1.1: Graph G = (V, E) with vertices V = {1, 2, 3, 4} and edges E = {(1, 2), (1, 3), (1, 4), (2, 4)}. The degrees of the vertices are d(1) = 3, d(2) = d(4) = 2 and d(3) = 1.

Figure 1 . 2 :

 12 Figure1.2: Distributions of the eigenvalues associated to matrices with order n = {10, 50, 100, 150, 500, 1000}, with iid Gaussian elements of mean 0 and standard deviation 2 in the diagonal, 1 otherwise. The red curves are related to Eq (1.27). According to Eq (1.27), the semi-circle law is an asymptotical behavior.

Figure 1 . 3 :

 13 Figure 1.3: Distributions of the eigenvalues associated to covariance matrices with order n = {10, 50, 100, 150, 500, 1000}, m = n, with iid Gaussian elements of mean 0 and standard deviation 1. The red curves are related to Eq (1.31). According to Eq (1.31), the Wishart distribution is an asymptotical behavior.

Figure 1 . 4 :

 14 Figure 1.4: Matrices with order n = {150, 500, 1000}. Each element is a realization of a Gaussian random variable with mean 10 and standard deviation 2 in the diagonal, 1 otherwise. The value of Λ 1 /Λ 2 are ∼ 60 if n = 150, ∼ 110 if n = 500, and ∼ 160 if n = 1000. When n increases the difference between Λ 1 [M ] and Λ 2 [M ] increases.

Figure 2 .

 2 Figure 2.2 plots another error number ||R|| versus αβ, for the four different types of matrix shown in the second column of Figure 2.1.The green dots refer to a realized ER graphs (with values of p from 0.3 to 0.8), the dark (light) blue dots refer to realized CB (BP) graphs (with k from 1 to n/2 -1), and the red dots refer to a realized SF graph (with γ from 1 to 8). As Figure2.2 shows ||R|| is low (high) for ER (SF) graphs for low values of αβ, and ||R|| increases when αβ increases. Whatever the values of p and γ, we always find a higher ||R|| for a SF than for an ER graph. However, the behavior of the CB graph is different: for small values of k (hence of αβ), the error is greater than the ones for ER graphs. Indeed, the size of both diagonal squares is almost the same. Since Λ 1 [M ] corresponds to the first square, and Λ 2 [M ] corresponds to the second one, the both eigenvalues are close. Hence for the same values of αβ, m (see Eq (2.27)) is higher for a CB graph than for an ER graph. Finally, when αβ increases, k also increases as indicated by the arrows in Figure2.2, and the size of the first diagonal square also increases. This means that when k increases, the checkerboard comes close to an ER graph, and Λ 2 [M ] decreases to 0 since the size of the second square decreases as well. However, the behavior of a BP graph is different from the CB graph: when k increases, αβ increases and ||R|| increases as well. Indeed, the non-diagonal upper block has increasing height, and decreasing width. Moreover, differently from a CB graph, |Λ 2 [M ]|/Λ 1 [M ] increases when αβ increases and then ||R|| increases. Bearing all above in mind, ER and SF graphs play an antogonist role, whereas CB and BP graphs play intermediary roles. In the right panels of Figure2.3, we show the behavior of the assortativity A D (see Eq (1.34)). For the CB and BP graphs, perfect networks are represented by the highest and lowest dots respectively. For each of these dots, the network has been 'rewired' preserving the degrees 100 times: one degree-preserving rewire as the following procedure. Randomly take two linked pairs of nodes, say (a,b) and (c,d). Break the links, and randomly associate a to c or d. Then associate b to d or c. This has the effect to decrease the absolute value of the assortativity, vertically moving on the up right plot of figure 2.2. However, the assortativity is not the only graph parameter which is related to αβ and ||R||. Looking at the CB graph, we see that the ability of the graph to be divided into modules changes the value of ||R||. When then to focus on the modularity. The modularity Q(G) of a graph G is the strength of a division of G into modules, and defined by the following equation:

Figure 2 . 1 :

 21 Figure 2.1: Four different simulated 200 × 200-adjacency matrices. First row: the ER graph. Starting with n unconnected distinct nodes, each pair of nodes is connected with probability p ∈ [0, 1]. The 1 are uniformly settled in the adjacency matrix, and the empirical histogram of the degrees approaches a Gaussian. The distance between the degree vector and the principal eigenvector is low. Second row: the SF graph. With an initial connected network of n 0 n nodes, n -n 0 nodes are added one at a time, and each of these new nodes are connected to n 00n 0 initial node(s) with a probability proportional to the initial degree of the n 00 nodes. The 'concentration' of 1's decreases when we go down right in the adjacency matrix. The distribution of the degrees of the nodes follows a power-law γ (the loglog histogram is also shown), and the distance between the degree vector and the principal eigenvector is larger than that of an ER graph. Third row: the CB graph. It is a graph splitted into two families. Each node is mainly connected with nodes of the same family with a probability p as in the ER graph's case, with rare connexion with the nodes of the other family. We then rewire 100 times preserving the degrees. The adjacency matrix is a chessboard, and the degree distribution is as two Gaussians. The distance between the degree vector and the principal eigenvector increases when the difference of the size of both families increases. Fourth row: The BP graph. It is a graph also splitted into two families. Each node of each family is exclusively linked with all the members of the other family. We then rewire 100 times preserving the degrees. The histogram is far from being Gaussian and the error is large.

Figure 2 . 2 :Figure 2 . 3 :Figure 2 . 4 :

 222324 Figure 2.2: Plot of the error ||R|| = nE 2 ( Ẽ) versus αβ, for the four graphs of Figure 2.1. For ER and SF graph types, the error ||R|| increases when αβ increases. In the CB graph, the error increases when αβ increases; the error has a maximum value, then decreases when αβ increases. Specifically when αβ increases from 1, |Λ 2 |/Λ 1 ∼ 1, and the error ||R|| increases to a maximum value (about 0.375). When αβ overcomes the point of the maximal value of ||R||, |Λ 2 |/Λ 1 decreases, since the first block size increases. The quantity |Λ 2 |/Λ 1 decreases faster than αβ inscreases: the error finally decreases towards 0. Contrary to the CB graph, αβ and |Λ 2 |/Λ 1 both increase in a BP graph, and ||R|| increases as well.

Figure 3 . 1 :

 31 Figure 3.1: 8000 × 8000 HiC matrix from human Embryonic Stem Cell (hESC), where the intra-chromosomal interactions have been removed (all the elements are null) (white diagonal blocks). This matrix is a checkerboard, and is then the object of clustering data analysis, leading to the identification of the nuclear compartments.

Figure 3 .

 3 Figure 3.4 shows an example of a checkerboarded matrix. The DNA contained in A highly interacts with itself, whereas DNA in B lowly does. An ideal case of HiC matrix would be that there is no any 1 in the AB and BA regions (see Figure3.5). Let M be the matrix of Figure3.5. We would like to compute its determinant, with a = a -λ and b = b -λ (which is the characteristic polynomial of a matrix with a = a and b = b). Let v 1 (resp. v 2 ) be the number of non-zero elements in any row of type A (resp. B), except the diagonal. v 1 and v 2 can be interpreted as the number of contact(s) that any part of the DNA in A (resp. B) has with itself in the same compartment. We immediately have

Figure 3 . 2 :

 32 Figure 3.2: Example of two nuclear compartments.The euchromatin compartment is the territory in the nucleus such that the DNA interaction with itself is low, and the heterochromatin compartment is the territory in the nucleus such that the DNA interaction with itself is high. © pixfocus.com

Figure 3 . 3 :

 33 Figure 3.3: Sum of each row of the 8000 × 8000 HiC matrix, from hESC, with respect to the principal eigenvector. The linear regression coefficient is 0.97

. 2 )

 2 Remark 3.2.1 Another method is to see that the characteristic polynomial does not change by permuting the rows and columns of M . In particular, it is possible to find a permutation P which transforms the matrix M into a blocked diagonal matrix P M t P , composed of two blocks. The first block is composed of elements a and a , the second of b and b . The whole determinant is the product of the determinant of each block. The first gives(a + v 1 a)(a -a) v 1 while the second gives (b + v 2 b)(b -b) v 2 ,which is exactly Eq (3.2). Moreover, thanks to the spectral theorem, it

Figure 3 . 5 :

 35 Figure 3.5: The considered ideal case of HiC matrix, of dimension n × n. The size of each rectangle is different. The elements in blue rectangles are all a ∈ R, in green ones b ∈ R, the diagonal in blue rectangles (red lines) are composed of a ∈ R, and in green ones (red dashed lines) b ∈ R. White rectangles are strictly composed of null elements. We are going to compute the eigenvalues of such a matrix.

. 3 )

 3 Setting now a = a -λ and b = b -λ, and a = b = 1 (adjacency matrix), and using (3.1), one has det(M -λId)

Figures 3 .

 3 Figures 3.6 and 3.7 plots some simulations of different HiC matrices. We distinguish two types of HiC matrix.

Figure 3 . 6 :

 36 Figure 3.6: HiC simulations: the matrix is a perturbed checkerboard, with some 1's in the AB and BA squares. The two highest eigenvalues are close, and the others tend to 0. Therefore V 1 is close to M e. However, there is no clear trend in the eigenvectors.

Figure 3 . 7 :

 37 Figure3.7: HiC simulations: the matrix is a perturbed checkerboard, with a very low number of 1's. The highest eigenvalue is far from the others, and V 1 is still close to M e. In this case, V 2 and also V 3 show interesting trends, and might be analyzed in deeper details to find the checker pattern.

Figure 3 . 8 :

 38 Figure 3.8: Finding the compartments from the HiC matrix : first enrich the data by a composition of type 1 and 2 (121 was the most revealing for the HiC data). The pattern appears inM 2 = Λ 2 t V 2 V 2 , whereas it is not clear for M 1 since V 1 is close to M e, which is noisy. We find the optimal number of clusters k (here k = 2) thanks to the criteria AIC, BIC and HQ, by doing successive kmeans in the (V 2 ,V 3 ) space. k is then the number of compartment(s), and we obtain the location of each DNA fragment in one compartment.

Figure 3 . 9 :

 39 Figure 3.9: One example of finding the compartments (with a simulated ABC chequered matrix (ABCAB)). The aim is to find in which compartment is located each DNA fragment. The procedure to do it is explained in Figure(3.8). One can see here that the optimal number of clusters is 3, corresponding to the number of compartments. The location of each DNA fragment is then found.

( 3 .

 3 10A)).

( 3 .

 3 10B)). As expected, there are more intra-compartmental contacts than inter-compartmental contacts, but the contacts in X inf are not existing: the removed rows and columns have indeed no 1's. Moreover the compartment A (resp. C) strongly interacts with the compartment B (resp. D), and the compartment A interacts more with the compartment D than with the compartment C. The loci of the genome occupying the compartments are presented in the Figure (3.10C).

Figure ( 3 .Figure 3 . 10 :

 3310 Figure(3.12) depicts the proportions of some biological features (row) over the compartments A, B, C and D (column). What we call 'proportion over A' is the ratio of the number of '1' at the same genomic position for the biological feature and the locations of A in the genome, out of the number of '1' at the same genomic position for the biological feature and the locations

Figure 3 .

 3 Figure 3.11: Matrix m whose element m ij is the proportion of paths from the compartment i to the compartment j obtained by following the whole genome from the extremity 5' to the other extremity 3' fragment by fragment, with i, j ∈ {A, B, C, D, X inf }. Specifically, the element m ij is the ratio of the number of paths from the compartments i to j, out of the number of paths from the compartments i to A, i to B, and so forth.

Figure 3 . 12 :

 312 Figure 3.12: Distribution of some biological features on the compartments. A The timing replication and Lamina Associated Domains (or LADs, the lamina protein presented on the nucleus surface, related to inactive DNA) pattern show that A and B are inactive compartments, whereas C and D are active compartments. B Proportion of Transcription Factors (TFs). C Proportion of histone marks.

Figure 3 . 13 :

 313 Figure 3.13: Proportions of the chromatin marks in each compartment. The chromatin marks responsible for the genome expression are mostly located in the fragments from the compartments C and D. The heterochromatin mark is mostly located in the compartment A. Repetitive sequences are located in the compartment X inf .

Figure 3 . 14 :

 314 Figure 3.14: Equivalent Figure 3.10 for the human Fibroblast cell. A The informational criteria gives that there are 2 nuclear compartments in the human fibroblast, named a and b. B Number of contacts between the compartments out of the genomic size in the compartments. b strongly interacts with intself whereas a less interacts with intself. C Distribution of the five compartments on the genome. underneath, the camemberg proportions of each of the hESC compartiments becoming the Fibroblast compartments.

  6 µm (Fig 3.15 b). From the bead coordinates we computed the associated Gram and distance matrices (Fig 3.15 c,d) and a binary contact map (Fig 3.15 a, steps 1 and 2 and Fig. 1e) and the distance matrix was then obtained by applying step 6 to this contact map (Fig 3.15 f). The consecutive application of steps 5 and 4 (Fig 3.15 a) reconstructs the coordinates up to a global transformation (some rotation, dilation and possibly mirror symmetry). To quantitatively assess the original structure recovery, we compared in a scatter plot the actual (Fig. 1d) and reconstructed (Fig 3.15 g) distances and computed their Spearman rank correlation (Fig 3.15 h).

  17 a), the chromosome average structure at 150-kb resolution (Fig 3.17 b) and the average arrangement of autosomal chromosomes within nuclear space at 3-Mbp resolution (Fig 3.17 c). The genome connectivity and chromosome partitioning achieved by ShRec3D (Fig 3.17

Figure 3 . 15 :

 315 Figure 3.15: Application of ShRec3D to a simulated data set. (a) Algorithm flowchart; steps 1 to 6 are detailed in the text and Online Methods. (b) 3D structure of in silico yeast genome generated using polymer dynamics for a chain of N = 26,538 beads (each chromosome is shown in a different color). (c) Gram matrix. (d) Distance matrix (s.u.). (e) Contact map (binary map, threshold = 60 nm). (f) Distance matrix derived from contacts. (g) Distance matrix of the reconstructed structure (dimensionless). (h) Scatter plot of original and reconstructed distances; heat map colors indicate the local density of points; Spearman rank correlation coefficient R is indicated. Inset, reconstructed 3D structure. s.u., simulation unit.

Figure 3 . 16 :

 316 Figure 3.16: Quantitative assessment of ShRec3D performance and reliability. (a,b) Comparison of ShRec3D with BACH, ChromSDE and MDS applied to inverse-frequency distances for simulated data of increasing size N (number of beads) in terms of reconstruction accuracy (Spearman rank correlation between original and reconstructed distances) (a) and computation time (b). (c) Robustness to a controlled amount of randomly misplaced contacts mimicking experimental noise (semilog plot). (d) Comparison of average distances in a population of an increasing number k of simulated structures (up to k = 500 independent snapshots of a Langevin dynamics of structure in Fig. 1b coarse-grained to N = 480 points) and distances reconstructed from the corresponding average contact map.

Figure 3 .

 3 Figure 3.17: 3D multiscale visualization of human autosomal chromosomes from Hi-C data. (a-f) Experimental contact maps (a-c) and corresponding 3D reconstructions (d-f) at genomic resolutions, from the scale of restriction fragments in chromosome 1 (embryonic stem cells (hESCs) (a) to that of bins containing 50 (b) or 1,000 (c) restriction fragments covering the whole chromosome set. Color gradients in d-f indicate the position along the genome.

Figure 3 .

 3 Figure 3.18: 3D-structure of the chromosome 1 in the hESC nucleus, obtained from the Fast-Floyd algorithm applied to the hESC intra-chromosomal HiC matrix. The compartments (A: blue; B: light blue; C: green; D: yellow) are the one obtained from the hESC inter -chromosomal HiC matrix already studied in this chapter. The big red ball is the centromer (associated with the compartment X inf ), and the two black balls delimit the chromosome. The parts of the chromosomes associated to a specific compartment are isolated from the others associated with another distinct compartment.

Proposition 4 . 3 . 5

 435 Let k pg < k pm s m /µ m . There exists τ 0 > 0 given by

  20 × 10 6 , k pm = 0.6, s m = 0.005, µ m = 0.002, k mp = 0.01, k pg = 2.95, s nr = 0.08, k nn = 0.01, k np = 0.1, µ nr = 0.12, k nn = 0.01, k np = 0.1, µ n = 0.05.

Figure 4 . 1 :

 41 Figure 4.1: In the left panel the nullclines of the pathogen (purple line) and of the activated phagocytes (black line). Zooming of the regions where there are the steady states: E 1 = (5 × 10 5 , 1.6) (center panel) and E 0 = (0, 0) (right panel).

Figure 4 . 2 :

 42 Figure 4.2: Numerical solutions of the delayed model (4.9)-(4.10) in the case τ 1 = τ 2 = τ and with initial conditions near the steady state E 1 = (5 × 10 5 , 1.6). The first column (a), (b) and (c) shows the time evolution of P (blue line), the second column (b), (e) and (h) shows the time evolution of N * (green line), the third column (c), (f) and (i) shows the associated phase space diagram of the delayed N * /P system (red line).

Figure 4 . 3 :

 43 Figure 4.3: Numerical solutions of the delayed model (4.9)-(4.10) in the case τ 1 = 0 and τ 2 0, and with initial conditions near the steady state E 1 = (50 × 10 6 , 1.6). The first column (a) and (d) shows the time evolution of P (blue line), the second column (b) and (e) shows the time evolution of N * (green line), the third column (c) and (f) shows the associated phase space diagram of the delayed N * /P system (red line).

  k pg =2.95, τ=0.4 k pg =2.7, τ=0.4

Figure 4 . 4 :

 44 Figure 4.4: Phase space diagram of the delayed model (4.9)-(4.10) in the case τ 1 = τ 2 = τ = 0.4. For k pg = 2.7 the system reaches the steady state E 0 (left panel) and for k pg = 2.95 the system reaches the steady state E 1 (right panel). The pathogen nullclines are in purple line and the activated phagocytes nullclines are in black line.

Figure 5 . 1 :

 51 Figure 5.1: The time evolution of the density of V (left panel) and of the density of KFc (right panel) for α = 0.2 and F 2 = 0 (solid line), F 2 = 0.0005 (dashed line), F 2 = 0.01 (dot line).

Figure 5 . 2 :

 52 Figure 5.2: The time evolution of the density of ISc for α = 0.5 and F 2 = 0 (solid line), F 2 = 0.0005 (dashed line).

Figure 5 . 3 :

 53 Figure 5.3: Distribution function of V for α = 0.8 and for F 2 = 0 (left panel) and for F 2 = 0.0005 (right panel).

Figure 5 . 4 :

 54 Figure 5.4: Distribution function of KFc for α = 0.8 and for F 2 = 0 (left panel) and for F 2 = 0.0005 (right panel).

Figure 5 . 5 :

 55 Figure 5.5: The time evolution of the density of KFc (left panel) and of the global activation energy of the system (right panel) for α = 0.2 and F 3 = 0 (solid line), F 3 = 0.0005 (dashed line), F 3 = 0.01 (dot line).

Figure 5 . 6 :

 56 Figure 5.6: Distribution function of the KFc for α = 0.2 and for F 3 = 0 (left panel) and for F 3 = 0.01 (right panel).

Figure 5 . 7 :

 57 Figure 5.7: The time evolution of the density of ISc for α = 0.8 and F 3 = 0 (solid line) and for F 3 = 0.001 (dashed line).

Figure 5 . 8 :

 58 Figure 5.8: Distribution function of the Cc for α = 0.8 and for F 4 = 0 (left panel) and for F 4 = 0.0005 (right panel).

Figure 5 . 9 :

 59 Figure 5.9: The time evolution of the density of ISc (left panel) and of the global activation energy (right panel) for α = 0.8 and F 4 = 0 (solid line), F 4 = 0.0005 (dashed line).

  Du ω(u * , u * ) du * = Ωu * ω(u * , u * ) du * = 1.

  3, . . . , n} (i.e. e is the unique eigenvector associated with 1). From t eN = t e, one has

	Hence
		t eV i [N ] = 0, ∀i ∈ {2, 3, . . . , n}.	(2.28)
	Let R be the column vector whose k-component reads:
	R k =	min j∈{1,2,...,n}	N kj =	min j∈{1,2,...,n} M jk , k ∈ {1, 2, . . . , n}.	(2.29)

  5.2.1 Since f i (t, u) du denotes the number of cells, at time t, such that u ∈ [u, u + du] then f

i (t, u * ) du * f j (t, u * ) du * (5.5)

is a non linear product that refers to the interaction, at time t, between the number of candidate cells such that u * ∈ [u * , u * + du * ] and the number of field cells such that u * ∈ [u * , u * + du * ]

  Du×Du η e ir B ir (u * , ω * , u) f i (t, u * ) g ir (t, ω * ) du * dω

	m
	η e ir g
	r=1 Du

* -f i (t, u) ir (t, ω * ) dω * ,

(5.22)

[M ] < 1, ∀j ∈ {1,

2, . . . , n}, then

t eN V i [N ] = t eV i [N ] ⇐⇒ (1 -λ i [N ]) t eV i [N ] = 0, ∀i ∈ {2, 3, . . . , n}.

-

The system is assumed to be homogeneously distributed in space. According to a mean field approximation, the encounter rate is assumed to be constant for all interacting pairs; for simplicity, η ij = 1 for all i, j ∈ {1, 2, 3, 4, 5}.

-

The term A ij related to the transition probability density is assumed to be defined by a delta Dirac function (deterministic output m ij (u * , u * ) of a pair interaction) depending on the microscopic state of the interacting pairs:

where m ij will be defined in the sequel. Under the above assumptions, the modelling of each term of interaction is obtained as follow.

• Interactions with transition of activity. It is assumed that V, KFc, and Cc are the only subsystems subject to conservative interactions, therefore

, and J 4 [f ](t, u) are derived under the assumption that the non-self cells have a tendency to increase their microscopic state with a certain rate, regulated by the interactions with the Fc and V. The evolution toward higher level of activity is called heterogeneity-proliferation for KFc and V instead is called heterogeneity-progression for Cc. Moreover, it is assumed that the heterogeneity-proliferation of KFc is greater than the heterogeneity-proliferation of V. According to Eq.(5.12), we define

where α is a positive parameter, which defines the heterogeneity-proliferation rate, toward high states of activity, for KFc, while ε < 1 is a scale parameter, which takes into account the difference among the heterogeneity-proliferation and the heterogeneity-progression rates of the V and Cc with respect to the KFc. Straightforward calculations give the following conservative terms:

(5.14)

• Interactions with proliferation-destruction processes. Each subsystem may proliferate when it encounters another subsystem. Specifically:

P.1 Fc proliferate when encounter each other;

P.2 V proliferate when encounter Fc and ISc; P.3 KFc proliferate when encounter Fc and V;

P.4 Cc proliferate when encounter V;

P.5 ISc proliferate when encounter V, KFc, and Cc.

Moreover, it is assumed that the proliferation rate of KFc cell is greater than the proliferation rates of Fc, V, and Cc; the proliferation rate of the ISc, when encounter V and Cc, is greater than the proliferation rate of the ISc when encounter KFc. Accordingly, we define:

where β is a positive parameter corresponding to the proliferation rate of KFc; β i is the proliferation rate of the ISc, and ε is the scale parameter. The proliferation terms thus read:

(

The cells of each subsystem may be destroyed when they encounter cells of another subsystem. Moreover, ISc are able to destroy (with rate δ) V and Cc more efficiently than the KFc; the latter are destroyed by V less efficiently than Fc. The non-self cells with an high level of activity have the ability to inhibit or destroy the ISc (immune suppression or immune-subversion) but the destruction rate of the ISc by KFc is assumed to be less than the destruction rate by Cc and V. It is worth stressing out that we assume a destruction of the ISc by non-self entities, since inhibited immune cells do not play a relevant role in the competition and may be equivalently assumed as eliminated. Accordingly, we define

where δ i is a positive parameter which characterizes the destruction of the ISc by V and Cc. Thus the destructive terms read:

(5.18)

In particular

• Interactions with mutation. It is assumed that may occur cell interactions which result, due also to genetic mutations, to the birth of a new subsystem. Specifically:

T.1 Fc may mutate in KFc when interact with Fc;

T.2 Fc may mutate in KFc when interact with V;

T.3 KFc may mutate in Cc when interact with V;.

In particular, we assume that is more likely that Fc become KFc when they encounter V. Moreover, the microscopic state of the entities does not change during the transition. Accordingly, we define

where γ and λ are the transition (mutation) rates in the KFc subsystem and in the Cc subsystem, respectively. Thus the transition terms read:

(5.20)

Finally we assume that a constant external force field F i , for i = {2, 3, 4}, acts on the ith functional subsystem and mimics a therapy.

Bearing all above in mind, the model consists in an evolution equation for each distribution function f i , for i ∈ {1, 2, 3, 4, 5}. Thus yields:

(5.21)

The model is characterized by 11 parameters having the following biological meaning:

• α is the heterogeneity rate of KFc;

• β is the proliferation rate of KFc;

• β i is the proliferation rate of ISc;

• δ is the destruction rate of V and Cc by ISc;

• δ i is the destruction rate of ISc by V and Cc;

• γ is the mutation rate of Fc in KFc;

• λ is mutation rate of KFc in Cc;

• is the scale factor.

• F i , for i = {2, 3, 4}, is an external force field that acts on the ith functional subsystem and mimics a therapy.

It is worth noting that the α-parameter refers to mass conservative encounters, the βparameters refer to encounters which generate proliferative events, the δ-parameters refer to