Raouf El Cheikh

Nader Samuel Bernard

El Khatib

Nader El Khatib
email: nader.elkhatib@lau.edu.lb

A multiscale modelling approach for the regulation of the cell cycle by the circadian clock

We present a multiscale mathematical model for the regulation of the cell cycle by the circadian clock. Biologically, the model describes the proliferation of a population of heterogeneous cells connected to each other. The model consists of a high dimensional transport equation structured by molecular contents of the cell cycle-circadian clock coupled oscillator. We propose a computational method for resolution adapted from the concept of particle methods. We study the impact of molecular dynamics on cell proliferation and show an example where discordance of division rhythms between population and single cell levels is observed. This highlights the importance of multiscale modeling where such results cannot be inferred from considering solely one biological level.

The mammalian cell cycle and the circadian clock are two molecular processes that operate in a rhythmic manner. On one hand, the cell cycle is driven by the rhythmic activity of cyclin-dependent kinases, which dictate the time a cell must engage mitosis and the time it must divide giving birth to two daughter cells. On the other hand, the circadian clock is a network of transcriptional and translational feedback-loops that generate sustained oscillations of different mRNAs and proteins concentration with a period of approximately 24h. It turns out that several components of the circadianclock regulate various cyclin-dependent kinases at different stages of the cell cycle. This makes the circadian clock a key player in the temporal organization of the cell cycle and makes these two biological processes act as two tightly coupled oscillators.

Many computational models explored the interaction between the cell cycle and the circadian clock. One category of models consists of systems of ordinary differential equations that describe the interaction at the molecular level [START_REF] Gérard | Entrainment of the mammalian cell cycle by the circadian clock: Modeling two coupled cellular rhythms[END_REF][START_REF] Zamborszky | Computational analysis of mammalian cell division gated by a circadian clock: Quantized cell cycles and cell size[END_REF]. The second category consists of physiologically-structured partial differential equations that describe the effect of circadian regulation on the growth rate of cells [START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF][START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF][START_REF] Bernard | Why do cells cycle with a 24 hour period?[END_REF][START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation[END_REF]. In a recent study, we have combined these two approaches by coupling an age-structured PDE to a circadian clock-cell cycle molecular ODE system [START_REF] Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF]. Although this model could capture the influence of intracellular dynamics on the growth rate of cells, it lacked a proper multiscale description [START_REF] Bernard | How to build a multiscale model in biology[END_REF]. It could not take into consideration, for example, intracellular heterogeneity among cells.

Here, we present a multiscale formulation of our mathematical model for the regulation of the cell cycle by the circadian clock and a numerical method for solving it.

The multiscale formulation consists of a transport equation structured by the molecular contents of the coupled cell cycle-circadian clock oscillator. In its general form, the equation reads ∂ t ρ(x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = L[x, λ] ρ(x, t, λ) ,

with ρ representing the density of cells. In the coupled PDE/ODE model, the PDE (equation [START_REF] Do | Wee1 kinase as a target for cancer therapy[END_REF] in [START_REF] Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF]) was structured in age. By contrast, partial derivatives in equation (1.1) are structured by the molecular contents (denoted x) and the space where it is solved is d-dimensional with d the number of molecular components x i ; this is why we call it a molecular-structured equation. Since the molecular mechanism of the cell cycle-circadian clock interaction involves an abundant number of components, the model is high-dimensional in nature. Classical numerical methods such as finite

A C C E P T E D M

A N U S C R I P T volumes/differences are inappropriate for solving the transport equation. This is the main difficulty that makes molecular-structured models scarcely used in similar applications. We circumvent this difficulty by adapting a particle method to solve the main equation.

Particle methods have arisen as an alternative to classical numerical methods for solving high-dimensional problems. They are used in different applications, like the incompressible Euler equations in fluid mechanics [START_REF] Harlow | The Particle-in-Cell Method for Fluid Dynamics[END_REF][START_REF] Leonard | Vortex methods for flow simulation[END_REF][START_REF] Leonard | Computing three-dimensional incompressible flows with vortex elements[END_REF], the Vlasov equation in plasma physics [START_REF] Bossy | A stochastic particle method for the McKean-Vlasov and the Burgers equation[END_REF][START_REF] Glassey | The Cauchy Problem in Kinetic Theory[END_REF] and in turbulence models for reactive flows [START_REF] Pope | Lagrangian PDF methods for turbulent flows[END_REF][START_REF] Pope | PDF methods for turbulent reactive flows[END_REF]. The computational implementation of particle methods is conceptually simple. At a given time, the solution is represented by a large number of particles, each having its own properties, for example position and weight. These properties evolve in time according to a system of stochastic or ordinary differential equations. The classical solution to the PDE is obtained by reconstructing a smooth density from the particles distribution.

Theoretically, the numerical solution is a linear combination of Dirac masses

ρ(x, t) = N j=1 α j (t)δ(x -x j )
where α j is the weight of particle j, x j its position (state) and N is the total number of particles. To obtain a classical solution, one has to update the positions and weights of particles and then regularize the Dirac masses. The overriding strength of particle methods is that, for N fixed, the size of the system increases only linearly with the dimensionality of the space. This means that if we have a structured equation with dimension d, we have to solve d × N ODEs/SDEs. Since recovering a classical solution requires a regularization of the Dirac distribution, the performance of particle methods depends also on the quality of the regularization procedure. For generalities about particle methods, one can refer to [START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Raviart | An analysis of particle methods[END_REF][START_REF] Cottet | Vortex methods: theory and practice[END_REF].

In this work, we investigated the influence of molecular dynamics on cell proliferation. Therefor, the novelty of our approach emanate from considering a multiscale design that incorporates three linked biological levels. Hereafter some specificities of our model:

• Intracellular heterogeneity among cells: generally speaking, the particle method consists of solving N copies of an ODE/SDE system, each representing the dynamics of a given particle. In our model, the N copies are not the same. Since each particle is associated with a given cell, we considered several sources of variability. For example, we have multiplied the

A C C E P T E D M A N U S C R I P T
right-hand side of the ODE system describing cell cycle dynamics by a factor λ. This ensures that each cell has a distinct cell cycle period.

Another source of variability was introduced with cell division. Once a cell divides, its daughter is attributed a new coefficient λ given by λ N +1 = λ N + Dξ λ , with D a "diffusion" coefficient and ξ λ a normally distributed random number. This implies that the daughter and mother cells will have distinct cell cycle periods.

• Inter-cellular connection: we assumed that cells are connected with each other by making Per/Cry mRNA transcription dependent on the average concentration of Per/Cry among cells. Including connectivity between cells is important in certain tissues to maintain a robust synchronized activity.

• Coupling back from population to molecular level : we assumed that cell division is dependent on the MPF-WEE1 molecular activity and that the total number of cells has an impact on Per/Cry mRNA activity.

This implies a two way coupling between the molecular and population levels. We assumed also that the MPF activation coefficient decreases at an exponential rate proportional to the total number of cells. This ensures a limited growth, which is physiologically more realistic.

The paper is divided mainly into two parts. First, we introduce the model structure,

give details about its construction and how the transport equation is solved by adapting a particle method. Then, we illustrate biological properties that can be examined using our model; like the effect of intracellular dynamics on cell proliferation, the synchronization among cells, and how heterogeneity among cells affects growth rate.

2 Description of the model

General framework and main equation

We consider a large collection of cells with state (x, λ) the parameter state is fixed over the lifespan of the cell, but can vary among cells or when a cell divides. The population is described by its density ρ(x, t, λ) which evolves according to the following hyperbolic PDE

∈ Ω × Λ ⊂ R d × R p ,
∂ρ ∂t (x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = 1 2 ∆ λ • σ 2 (x, t, λ)R(x, λ)ρ(x, t, λ) + r(x, λ)ρ(x, t, λ) L[x, λ] ρ(x, t, λ) . (2.2)
We assumed that the cellular state dynamics x is governed by a deterministic system of ODEs, represented by the term u(x, t, λ, ψ), which depends on the cellular state (x,λ)

and on m population-level statistics ψ : Ω × R × Λ → R m where

ψ i = ρ, Φ i = Ω×Λ ρ(x, t, λ)Φ i (x, t, λ)dxdλ, i = 1, • • • , m (2.3) with Φ i : Ω×R×Λ → R, i = 1, • • • , m taken such that ψ i is finite. The operator L[x, λ]
describes the population relative growth rate. The terms σ, R and r are diffusion in parameter state, differentiation and growth rates of the population, respectively. For the sake of simplicity, we omitted the full description of the derivation of equation (2.2).

The reader could refer to Supplementary materials for more details. We adopted a more computational and biological description of the problem.

Intracellular dynamics

The evolution of cells in the state space depicts the cell cycle-the circadian clock coupling mechanism. Therefore, intracellular dynamics were described according to the following deterministic ODE system presented in [START_REF] Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF]:

dx dt = u(x, t, λ, ψ). (2.4) A C C E P T E D M A N U S C R I P T
where

u 1 = dx 1 dt = ν 1b (x 7 + c) k 1b (1 + ( x 3 k 1i ) p ) + x 7 + c -k 1d x 1 + k s ψ 1 , (2.5
)

u 2 = dx 2 dt = k 2b x q 1 -k 2d x 2 -k 2t x 2 + k 3t x 3 , (2.6
)

u 3 = dx 3 dt = k 2t x 2 -k 3t x 3 -k 3d x 3 , (2.7
)

u 4 = dx 4 dt = ν 4b x r 3 k r 4b + x r 3 -k 4d x 4 , (2.8
)

u 5 = dx 5 dt = k 5b x 4 -k 5d x 5 -k 5t x 5 + k 6t x 6 , (2.9
)

u 6 = dx 6 dt = k 5t x 5 -k 6t x 6 -k 6d x 6 + k 7a x 7 -k 6a x 6 , (2.10)
u 7 = dx 7 dt = k 6a x 6 -k 7a x 7 -k 7d x 7 , (2.11
)

u 8 = dx 8 dt = λ k lmpf +k 0mpf exp -ηψ 2 k n 1mpf k n 1mpf +x n 8 +sx n 10 (1 -x 8 ) -d wee1 x 9 x 8 , (2.12
)

u 9 = dx 9 dt = λ k actw k actw + d w1 (c w + Cx 7 ) + k actw k actw + d w1 -1 k inactw x n 8 x 9 k n 1wee1 + x n 8 -d w2 x 9 , (2.13)
u 10 = dx 10 dt = λ k act (x 8 -x 10 ) . (2.14)
Equations (2.5-2.11) describe the circadian oscillator and equations (2.12-2.14) the cell cycle. Circadian dynamical variables are x 1 , Per2 or Cry mRNA and proteins;

x 2 , PER2/CRY complex (cytoplasm); x 3 , PER2/CRY complex (nucleus); x 4 , Bmal1 mRNA; x 5 , BMAL1 cytoplasmic protein; x 6 , BMAL1 nuclear protein and x 7 , active BMAL1. Cell cycle dynamical variables are x 8 , active MPF; x 9 , active WEE1 and

x 10 , active MPF inhibitor. The coupling between these two oscillators is taken into consideration in the term Cx 7 of equation (2.13) through the regulation of WEE1 by BMAL1/CLOCK. The coefficient C is the coupling strength. Details about this model can be found in [START_REF] Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF]. The cellular dynamics is now part of a multiscale system, some changes were added to the original system to reflect that cells interact with each other.

These changes are detailed in the coming paragraphs.

Population synchronization

Including connectivity between cells is important in certain tissues to maintain a robust synchronized activity. It is known, for example, that neuronal clocks within the the suprachiasmatic nucleus SCN form a heterogeneous network that must synchronize to

A C C E P T E D M A N U S C R I P T
maintain time keeping activity. The coherent output of the SCN is established by intracellular signaling factors, such as vasointestinal polypeptide [START_REF] Enright | Temporal precision in circadian systems: A reliable neuronal clock from unreliable components?[END_REF]. A simple way to induce synchronization in our model is to make Per/Cry mRNA transcription depends on the average concentration of Per/Cry among cells. For that, we included in equation (2.5) the term k s ψ 1 := k s ρ, Φ 1 with k s a coefficient that represents the connectivity strength and Φ 1 = x 1 . Although peripheral tissues, as modeled here, would run outof-phase without the SCN, there are so many factors that can synchronize clocks that it is relevant to study the possibility for peripheral tissues to be connected.

Variability among cells

Variability among cells arises naturally from the difference in their molecular contents.

This was taken into consideration in our model through the initial states of cells, which can be chosen randomly. We added another source of variability through the parameter λ. We assumed that each cell has a distinct cell cycle period. This was modeled by multiplying equations (2.12-2.14) by a scaling factor λ. Cell division can also induce variability by assigning to each new born cell an intrinsic cell cycle period that is different from its mother cell. (This is accomplished by taking a non trivial distribution p(λ|z), see Supplementary materials equations (1.4), (1.5)).

Cell division

A central aspect of our model is that it takes into account cell division. Cell division timing is based on the antagonistic relationship between the mitosis promoting factor MPF and the protein WEE1. It is assumed that a cell enters mitosis once the activity of MPF surpasses that of WEE1 and then divides once MPF activity shuts down abruptly.

Based on this mechanism, we considered two division types, one deterministic and one stochastic. The deterministic division occurs exactly every time MPF activity rises above WEE1 activity and then shuts down. This was taken into consideration with a growth rate r(x, λ) := r d (x, λ) = δ (x∈Γ) . For the stochastic division, we considered that a cell divides with a certain probability t × r(x, λ) + o( t) for a small time step t. The division rate r is a function of x and λ that mimics the deterministic case.

For example, a switch-like function that takes small values on one side of Γ and large values on the other side could be used (see Supplementary materials for more details about the growth rates and Γ). Here, we used the Koshland-Goldbeter function [START_REF] Goldbeter | An amplified sensitivity arising from covalent modification in biological systems[END_REF] A

C C E P T E D M A N U S C R I P T
given by

K(y, z) = 2yJ i z -y + zJ a + yJ i + (z -y + zJ a + yJ i ) 2 -4yJ i (z -y) . (2.15)
This function generates a switching behavior [START_REF] Novak | Mathematical model of the cell division cycle of fission yeast[END_REF]. If the ratio y/z (y and z are dummy variables here) becomes larger than one, the function switches to the upper state and the transition occurs. J a and J i are two constants that determines the stiffness of the switch, if they converge to zero, the switch converges to a step function.

Limited growth

To make sure the growth is bounded (for physiological and computational reasons), we assumed that cell proliferation slows down when the total number of cells reaches a threshold value. We assumed that the activation coefficient of MPF decreases at an exponential rate proportional to the total number of cells. This was taken into consideration in equation (2.12) by multiplying MPF activation coefficient k 0mpf by exp -ηψ 2 where ψ 2 = ρ, Φ 2 with Φ 2 = 1 (ψ 2 = N (t)). When cell number is large enough, MPF activity cannot increase above that of WEE1, and cell division is blocked.

Method of resolution

The main difficulty for solving equation (2.2) comes from the high dimension of the dynamical state space (10 in our case). For that, we adapted a particle method that is more appropriate in high dimensions than classical methods such as finite differences/volumes.

The intuition behind the particle method is to start with a distribution of particles that approximates the initial condition and then follow the evolution of these particles in time in the dynamical and parameter states space. Mathematically, the particle solution is a discrete measure function, which means it is irregular. To obtain a solution in the usual classical sense at a given time T , one has to recover the classical solution with regularization techniques. The method is mainly divided into two steps

Step 1: approximation of the initial condition

Given an initial condition ρ 0 (x, λ) ∈ C 0 (R d × R), we take an initial set of particles x j with weights α j such that

ρ 0 h = N j α j δ[(x, λ) -(x j , λ j )] approximates ρ 0 . This should A C C E P T E D M A N U S C R I P T
be understood as an approximation in the sense of measures, which means that one looks for a test function

φ ∈ C 0 0 (R d × R) such that ρ 0 , φ = R d ρ 0 φdxdλ and ρ 0 h , φ = j α j φ(x j , λ j ).
This yields a typical numerical quadrature problem where the integral

R d ρ 0 φ is ap- proximated by N j
α j φ(x j , λ j ). For biological purposes, in our study we do consider directly an initial condition of the form (α j = 1)

ρ 0 h = N j δ[(x, λ) -(x j , λ j )].

Step 2: particle evolution in time

In the dynamical state space the particles positions X j (t) (not to confuse between X j and x j , our notations implies that X j (0) = x j ) can be computed at a given time t by solving the following ordinary differential equation system

d dt X j (t) = u(X j (t), t). (2.16)
To take into account the source term in equation (2.2), we examined the intracellular state of each cell (position of X in the state space) after each time step t and looked if the biological conditions are met for a cell to divide. This means that if for a given particle, MPF activity rises above that of WEE1 and then shuts down under a threshold level, we add a particle X N +1 with intracellular state similar to the mother cell or randomly chosen. If division occurs, a new parameter λ N +1 = λ N + Dξ λ is assigned to the particle of number N + 1, with D a diffusion coefficient and ξ λ a normally distributed random number. Hence the measure solution at a given time T , is given by

ρ h (x, T ) = Ñ j δ[(x, λ) -(X j (T ), λ j (T ))].
(2.17)

with Ñ the new number particles.

Commonly, with particle method one needs to obtain a solution in the classical sense, which means that the particle solution (2.17) should be

A C C E P T E D M A N U S C R I P T

regularized. This is usually accomplished by taking the convolution product of ρ h with a smoothing kernel ζ . In our case, we did not need to obtain a smooth solution, since there is no explicit dependence of the quantities of interest on the density of cells ρ. We focused on studying statistic-like properties of the distribution of cells. Hence, we did not perform any regularization procedure on the particle solution. For more information about particle methods the reader could refer to [START_REF] Pope | PDF methods for turbulent reactive flows[END_REF][START_REF] Raviart | An analysis of particle methods[END_REF].

3 Results

Regulation of the cell cycle by the circadian clock

We solved equation (2.2) with L = 0, i.e. without division. We assumed that each cell has an intrinsic cell cycle period included in the interval [20 h, 28 h] . This was done by assigning to each particle a different coefficient λ. We let all cells have the same initial molecular state and considered two cases; one without coupling to the circadian clock (C = 0), and one with coupling to the circadian clock (C = 1.2). We plotted a projection of the solution in the subspaces (x 7 , x 9 ) and (x 7 , x 8 ), which are the subspaces of molecular concentrations (BMAL1/CLOCK, WEE1) and (BMAL1/CLOCK, MPF).

We obtained a random distribution of particles for the non-coupled case and a limit cycle distribution for the coupled case (Figure 1). This means that the circadian clock is forcing all cells to oscillate with the same period. This is in agreement with results obtained earlier [START_REF] Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF]; it was shown that a population composed of cells with cell cycle periods between 20 h and 28 h are brought to oscillate with a unique cell cycle period of 24 h for large values of C.

Cell division

We took four cells initially, two of them with λ = 230, and two with λ = 180 (intrinsic cell cycle lengths 10 h and 12.7 h respectively). These values were taken randomly.

We solved equation (2.2) with a source term and considered three cases, one with r = r d (deterministic growth) and the two other ones with r = r s (stochastic growth). For the stochastic growth, we used a stiff and a non-stiff Koshland-Goldbeter function. First, we fixed the value of η to 0, which means that MPF activity does not depend on the total cells number. We assumed that newborn cells retain the cell cycle length of their mother cells at the division time, hence p(λ|z) = δ(λz) which leads the growth term to ne L = r(x, λ)ρ. This ensures that MPF cycle does not change along birth. Our

A C C E P T E D M A N U S C R I P T

simulations showed that the total number of cells after four days was equal to 1152 for the first case, 900 for the second one with a stiff Koshland-Goldbeter function and 480 with a non stiff one. This is natural and is justified by the fact that the division rate r d is a Dirac-type division which means that division depends deterministically on MPF activity. Looking at MPF activity, we see that it peaks 9 times for λ = 230, and 6 times for λ = 180 every 4 days (Figure 2A,B). If division occurs exactly once MPF activity accomplishes a normal cycle, the total number of cells should be 2 × 2 9 + 2 × 2 6 = 1152, which is the result obtained with r = r d (Figure 2C). The division rate r s introduces stochasticity in the decision for division. In this case, division does not depend only on MPF-WEE1 activity, but also on the probability r s × t at each time step t Figure (2C).

Second, we assumed that proliferation depends on the total number of cells, and set η = 2 × 10 -3 (increasing the value of η will make the factor exp -ηψ 2 less than one and hence decreases the value of MPF activation coefficient k 0mpf ). We took initially 100 cells with different intrinsic cell cycle periods and followed the total number of cells over 15 days. We remarked that division stopped after 8 days, and the growth curve reached a steady state (Figure 2D). Increasing the value of η decreases the activity of MPF with increasing number of cells. This means that at a certain time, MPF activity will decrease below a threshold that does not allow the cell to start mitosis.

Synchronization

To study the synchronization of cells, we compared the molecular concentrations of a cell chosen randomly with the average molecular concentration of all cells. We let the rate of production of Per2/Cry mRNA of each cell to be dependent on the average value for Per2/Cry mRNA. This was done by taking k s = 0.05 in equation (2.5). We took a population of 100 cells with autonomous cell cycle periods randomly distributed between 20 and 28 hours. Initial molecular concentrations were chosen randomly between 0 and 1 and each cell has a distinct initial molecular state. We did not consider division hence L = 0, but we considered coupling between the cell cycle and the circadian clock (C = 1.2). We followed the evolution of three components, Per2/Cry mRNA, BMAL1/CLOCK and WEE1, over 20 days. Our simulations showed that for k s = 0, the average molecular concentration had small oscillations which are asynchronous with those of the random cell concentration (Figure 3 A,B,C). Whereas for k s = 0.05, we obtained that the average molecular concentration tends to coincide with that of the random cell (Figure 3 D,E,F). This indicates that all cells are oscillating in a similar manner (same period and phase), and indicates synchronization of all popu-

A C C E P T E D M A N U S C R I P T

A C C E P T E D M

A N U S C R I P T lation cells. In the case k s = 0, even though cells are coupled in the same manner to their circadian clock, cells keep oscillating in an asynchronous manner. The circadian clock regulates all cells to oscillate in a similar period which is 24 h in this case but due to the randomness in the initial molecular concentrations, each cell oscillates with a different phase.

Heterogeneity of cells and growth rate

We studied the growth rate of a cell population where each cell had a different cell cycle period. We took 100 cells initially and conferred to each of them a coefficient λ chosen randomly between 128 and 193. These coefficients scale the intracellular cell cycle system so that periods range randomly from 12 h to 18 h. The circadian control strength value C was fixed to 1.6. We did not consider connection between cells (k s = 0) and we did not consider dependence of the intracellular dynamics on any population level statistics. Simulations were done over 20 days and with a net death rate to mainly reduce the computation time. Simulations showed a growth rate with two daily peaks, suggesting that there were two division rounds every day (Figure 4A). To explain this bimodal behavior for the growth rate, we looked at the average activity of WEE1 and MPF. Simulations showed that MPF activity overcomes that of WEE1 once a day (Figure 4B). Based on average MPF activity alone, division should occur only once a day, and cannot explain the two daily peaks. We then followed a subpopulation of cells and examined at which time each cell divided. We identified three regimes for division: a single division per day, two divisions per day and three divisions every two days (Figure 4C,D). These regimes can be explained by the entrainment results obtained previously (for autonomous cell cycle periods between 12 and 18 h) showing large phase-locking regions entrainment with ratios 1:1, 2:1 and 3:2 [START_REF] Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF]. The double peaks for the growth rate can be then explained by the fact that all cells are dividing at least once a day and some of them will undergo a second round of division (Figure 4D). All of these works emphasize on the importance of taking into account different levels of a biological system by showing results that cannot be inferred by considering solely one level of functioning. We exposed a simulation showing discordance between rhythms observed on the population and signal cell levels. This is important since it underlies the fact that biological function markers which most of the times rely on average-like information can mislead the interpretation in some cases. Our modeling approach differs from the above cited works which can be separated into two categories. One consisted of low dimension structured partial differential equations solved by classical numerical techniques [START_REF] Ribba | A multiscale mathematical model of cancer, and its use in analysing irradiation therapies[END_REF][START_REF] Brikci | Analysis of a molecular structured population model with possible polynomial growth for the cell division cycle[END_REF]; and the other one consisted of agent-based modeling approaches [START_REF] Bratsun | Multiscale modeling of tumor growth induced by circadian rhythm disruption in epithelial tissue[END_REF][START_REF] Schlüter | Multi-scale modelling of the dynamics of cell colonies: insights into cell-adhesion forces and cancer invasion from in silico simulations[END_REF][START_REF] Shokhirev | A multi-scale approach reveals that nf-b crel enforces a b-cell decision to divide[END_REF]. Our approach consists of considering a master partial differential equation in high dimension structured by molecular contents for which we propose a computational solution based on a particle method. This allowed us to structure the partial differential equation by a large number of molecular contents, providing a realistic representation of the coupled cell-cycle circadian clock oscillator.

Discussion and conclusion

By doing so, we were able to describe both population and single cell functioning in one equation. Globally, the model takes into consideration the molecular dynamics inside a cell, connection between heterogeneous cells and the proliferation of cells.

Our results can shed light on the superiority of frequent administration of chemotherapies over the standard maximum tolerated dose (MTD). Scheduling the MTD into

A C C E P T E D M

A N U S C R I P T multiple administrations per cycle proved to be more effective and less toxic [START_REF] Hénin | Revisiting dosing regimen using pk/pd modeling: the model1 phase i/ii trial of docetaxel plus epirubicin in metastatic breast cancer patients[END_REF][START_REF] Pasquier | Metronomic chemotherapy: new rationale for new directions[END_REF][START_REF] Slevin | A randomized trial to evaluate the effect of schedule on the activity of etoposide in small-cell lung cancer[END_REF].

We showed in section 3.4 that an apparent population growth rhythm may hide several subpopulation growth rhythms. Biological markers, as the Ki-67 index, are mostly indicative of the growth status at the population level. Therefore, one is tempted to consider that the administration of an anticancer drug during the time of the largest growth peak is the most effective. However, our simulations showed that this largest growth peak is obtained from the crossing of different subpopulation rhythms. There-

  S C R I P T Highlights • A multiscale model for the cell cycle-circadian clock coupling is developed • The model is based on a transport PDE structured by molecular contents • A particle-based method is used for resolution • Impacts of inter and intracellular dynamics on cell proliferation are studied • Discordance of division rhythms between population and single cell levels is observed A C C E P T E D M A N U S C R I P T A multiscale modelling approach for the regulation of the cell cycle by the circadian clock A C C E P T E D M A N U S C R I P T 1 Introduction

  where x is a cellular dynamical state in an open subset Ω of dimension d and λ a cellular parameter state in an open subset Λ of dimension p. In our case, x = x 1 , • • • , x d with d = 10 represents the set of proteins/mRNAs concentrations and λ ∈ R (p = 1) represents the intrinsic cell cycle period of each cell. The distinction between the dynamical and the parameter state is that the dynamical state changes during the lifespan of a cell, while A C C E P T E D M A N U S C R I P T

Figure 1 :

 1 Figure 1: Solution of equation (2.2) without cell division. (A,B) Projection on the subspace (x 7 , x 9 ); (C,D) projection on the subspace (x 7 , x 8 ). (A,C) Without coupling to the circadian clock, C = 0. (B,D) With coupling to the circadian clock C = 1.2.

Figure 2 :

 2 Figure 2: (A) MPF-WEE1 antagonistic activity for λ = 180: division occurs if MPF activity rises above that of WEE1 and then shuts down. (B) Impact of λ on MPF activity. (C) Total number of cells with different source terms with no dependence of MPF activity on the total number of cells (η = 0). Initial number of cells = 4, two cells with λ = 230 and two with λ = 180. (D) Impact of the dependence of MPF activity on the total number of cells (η = 2 × 10 -3 ).

A

  We presented in this work a multiscale mathematical model for the regulation of the cell cycle by the circadian clock and its influence on cell growth. The model consisted of a master transport partial differential equation of high dimension structured by the molecular contents of the coupled cell cycle-circadian clock oscillator. We proposed a

Figure 3 :Figure 4 :

 34 Figure 3: Comparison between population average molecular concentrations randomly chosen and single cell molecular concentrations. (A,B,C) No synchronization between cells: k s = 0; (E,F,G) Synchronized activity: k s = 0.05.

  fore, administrating a drug at the same time every treatment cycle may hit only one subpopulation of cells, leaving another subpopulation without treatment for a long time. This possibly leads to the emergence of resistant subpopulations and make the treatment ineffective.The aim of this work was to provide a computational multiscale framework for the regulation of the cell cycle by the circadian clock. The use of multiscale computational tools is of growing interest in the area of drug development, especially the area of systems pharmacology. A major challenge for systems pharmacology is the development of mechanistic tools to understand how regulatory networks control variability in drug response at the organismal level[34]. Our model was constructed in the hope to overcome such challenges. It had the specificity of describing the connection between the cell cycle and the circadian clock through the regulation of the protein kinase WEE1 by the complex BMAL1-CLOCK. Irregular activity of the WEE1 kinase has been linked to ovarian and NSCLC cancers. Inhibiting this activity demonstrated interference with DNA damage response within tumor cells, leading to their death. Drugs targeting WEE1 are currently under investigation[START_REF] Do | Wee1 kinase as a target for cancer therapy[END_REF]. In such specific applications, our model can be of great utility. According to Zhao and colleagues, experimental settings in the process of drug development are sometimes not sufficient due to the multitude of interactions between mammalian network proteins. Actions on the desired targets may lead to adverse events occurring at another level due to the wide ramification of molecular interactions. The abundance of such interactions makes it hard to design drugs affecting only the desired targets with controlled therapeutic effects[34]. Therefore it is important to propose modeling frameworks that allow in silico experimentation in the presence of uncertain or incomplete data. In this manuscript, we have provided such an example, and showed why incomplete population level data may lead to a suboptimal treatment.Multiscale modeling in systems medicine is growing in importance. Wolkenhauer and colleagues recently proposed different recommendations to improve the integration of multiscale models into clinical research[START_REF] Wolkenhauer | Enabling multiscale modeling in systems medicine[END_REF]. One of their medium-term recommen-A C C E P T E D MA N U S C R I P T

A C C E P T E D M A N U S C R I P T

dations is the rising need of developing computational tools and algorithms for efficient multiscale simulations. We hope the modeling approach we propose in this study represents a step forward towards this direction.

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T

[34] S. Zhao and R. Iyengar. Systems pharmacology: network analysis to identify multiscale mechanisms of drug action. Annu. Rev. Pharmacol. Toxicol., 52:505-521, 2012.


Bastien Durand

Franck Delvare

Patrice Bailly

Didier Picart

Identification of the friction under high pressure between an aggregate material and steel: experimental and modelling aspects

The interface behaviour between steel and a quasi-brittle aggregate material is characterised up to normal pressures of the magnitude of 100MPa. This article presents a new test enabling the behaviour of the interface to be studied whilst retaining the sample's integrity. The experimental configuration having been retained consists in sliding a cylindrical sample of the material inside a steel tube, said tube acting both as a sliding surface and containment ring. The sample is pushed on one side and faces a spring on the other. The axial compression generates the interface pressure by Poisson effect. This originality of this assembly lies in the simultaneous application of normal pressure to the interface and of its relative motion. The assembly is placed in a quasi-static testing machine. The analysis is made by means of an analytical modelling of the test. This method enables the identification of the initial contact conditions (tightening of the sample), the friction coefficient and its dependence on the pressure. Numerical simulations of the test using a finite element method enables the analytical approach to be validated as well as the set of parameters identified depending on the normal pressure regimes.

Introduction

Certain dynamic tests enable the reactivity of compressed explosives (PBX) subjected to shocks to be tested. These are namely: the drop-weight test [START_REF] Field | Ignition mechanisms of explosives during mechanical deformations[END_REF], the Steven-test [START_REF] Gruau | Ignition of a confined high explosive under low velocity impact[END_REF], [START_REF] Vandersall | Experimental and modeling studies of crush, puncture, and perforation scenarios in the Steven impact test[END_REF] and the Taylor shock test. The place and time of the priming of the chemical reaction particularly depend on the conditions of contact with the wall [START_REF] Dickson | Frictional Heating and Ignition of Energetic Materials[END_REF], [START_REF] Hoffman | Aspect of the tribology of the plastic bonded explosive LX-04[END_REF], [START_REF] Picart | Ignition of a HMX-based PBX submitted to impact: strain localisation and boundary condition[END_REF]. A numerical simulation of one variant of the Taylor test forms a simple illustration of the influence of the friction conditions upon ignition. This test consists in projecting a sample of explosive (cylindrical sample with a diameter and a height of 18mm) against a steel wall assumed to be pressure-resistant. The thermomechanical characteristics of PBX are given in [START_REF] Picart | Non-shock ignition of a HMX-based high explosive: thermomechanical numerical study[END_REF]. Figure 1 shows the influence of the friction coefficient f on the maximal temperature reached by the sample during the impact. A very significant difference of around 200K is obtained between the two cases.

During an impact at these velocities, the normal contact pressures and the sliding velocities between steel and explosive are respectively of around 100MPa and 10 to 100m/s.

Our aim is to recreate these velocity and pressure conditions, but in a configuration enabling the friction to be measured. Certain test benches partially fulfil this aim: tribometer with explosively-generated friction [START_REF] Kim | Nanostructures generated by explosively driven friction: Experiments and molecular dynamics simulations[END_REF], target-projectile assembly with oblique impact [START_REF] Rajagopalan | Novel experimental techniques for investigating time resolved high speed friction[END_REF], torsional Hopkinson bars [START_REF] Huang | Dynamic tribological response of SiC fracture surfaces[END_REF], [START_REF] Huang | Dynamic Friction of SiC Surfaces: A Torsional Kolsky Bar Tribometer Study[END_REF], [START_REF] Rajagopalan | Novel experimental techniques for investigating time resolved high speed friction[END_REF], [START_REF] Rajagopalan | An experimental method to study high speed sliding characteristics during forward and reverse slip[END_REF], dynamometrical ring with parallelepipedic sample launched by a gas cannon or by a hydraulic machine [START_REF] Philippon | An experimental study of friction at high sliding velocities[END_REF] and possibly pin on rotating disk [START_REF] Dickson | Frictional Heating and Ignition of Energetic Materials[END_REF].

For the sake of safety, the experiments are developed using an inert material, denoted I1. This material is namely a fine compressed powder (mixture of barium meal BaSO 4 -29.3%wt -, of melanime -65.2%wt -and a binder -5.5%wt -made of an epoxy resin). Its mechanical behaviour closely resembles that of an explosive. This behaviour has previously been studied by carrying out triaxial compression tests [START_REF] Bailly | Dynamic behavior of an aggregate material at simultaneous high pressure and strain rate: SHPB triaxial tests[END_REF]. Under compressive loading with high hydrostatic pressure (which is different from a purely hydrostatic loading), the material is able to flow when its plasticity threshold has been attained (here the maximal constraints obtained using triaxial tests are assimilated to a plasticity threshold in order to simplify the behaviour model). The plasticity threshold thus identified is of the Drucker-Prager type [START_REF] Bailly | Dynamic behavior of an aggregate material at simultaneous high pressure and strain rate: SHPB triaxial tests[END_REF].

That is σ, P and σ mises respectively the stress tensor, the hydrostatic pressure and the Von Mises equivalent stress:

( where tr is the trace operator and I the identity tensor. The criterion is thus given by the following formula:

(2)

C P mises < -α σ
where the stress C depends on the strain rate. It equals 25MPa at low strain rates and exceeds 80MPa at those strain rates attained during compression tests on Hopkinson bars (typically 100s -1 ). The coefficient α is constant and equal to 0.64. In the event of uniaxial compressive loading, the criterion corresponds to a maximal acceptable stress of around 30MPa for low strain rates and around 100MPa for strain rates reached on Hopkinson bars [START_REF] Bailly | Dynamic behavior of an aggregate material at simultaneous high pressure and strain rate: SHPB triaxial tests[END_REF].

One way of subjecting the interface to high pressures without fracturing the sample is to confine it. For this, a cylindrical sample of I1 is encased in a quasi-rigid tube. This technique is employed to perform compression tests with quasi-uniaxial strain states [START_REF] Bailly | Dynamic behavior of an aggregate material at simultaneous high pressure and strain rate: SHPB triaxial tests[END_REF], [START_REF] Forquin | An experimental method of measuring the confined compression strenght of geomaterials[END_REF], [START_REF] Forquin | Influence of free water on the quasi-static and dynamic of strength of concrete in confined compression tests[END_REF]. Friction between the tube interface and the sample is here perceived as a drawback and has been studied in [START_REF] Azhdar | Determination of dynamic and sliding friction, and observation of stick-slip phenomenon on compacted polymers during high velocity compaction[END_REF], [START_REF] Burlion | Experimental analysis of compaction of concrete and mortar[END_REF], [START_REF] Tien | Wall Friction measurement and compaction characteristics of bentonite powders[END_REF]. Since the mechanical behaviour of the tube is known, the strain measurements made by means of gauges glued to the external face enable the stresses on the interface to be measured. The idea retained is to slide the cylindrical sample into the tube, which acts as both a sliding surface and a confining vessel (given the Poisson effect of the sample). This principle has already been tested with a hollow sample, through whose length a tightening screw [START_REF] Durand | Caractérisation du comportement à l'interface entre un matériau agrégataire et un acier sous forte pression et à grande vitesse. 20 ième Congrès Français de Mécanique[END_REF] passes. Despite the advantages of this principle, this first device does not allow the sample to be subjected to high pressures. The device being studied in this article allows us to overcome this drawback and to reach pressures of 100MPa at the interface. Our study here is limited to the cases of low relative velocities between the bearing faces.

The detailed description of the experimental configuration and the limits linked to the internal sliding are presented in section 2. A model enabling the identification of the friction parameters based on the measurements is described in section 3. The experimental results on the test examples using an inert material are then presented and analysed in section 4. Finally, section 5 focuses on validating the analysis method.

The experimental layout

The experimental device

The device proposed has been designed to reach high pressure levels. Before the test, the cylindrical sample is confined in a steel tube (c.f. Figure 2). It is pushed from top to bottom by a mechanical machine. The movement of its base is countered by a spring assembly (of the Belleville type). The stiffness of the spring allows the value of the counterpressure to be controlled. The originality of this assembly lies in that the normal pressure (by Poisson effect on the sample) and the sliding are imposed simultaneously. This allows the full pressure range to be covered, with the pressure gradually increasing along with the displacement. The inner wall of the steel tube was reamed and the sample was turned on a sliding lathe. Both have a weak surface roughness.

Figure 3 is a photograph of the device. For practical reasons, the confinement tube and the tube enclosing the spring may be screwed together. The parts belonging to such fastening system are neither shown nor modelled in the simulations by the finite element method mentioned in section 5.

The values of forces F m and F r are measured during the tests (c.f. Figure 4). F m is the force exerted by the machine on the sample. It is given directly by the machine's force sensor.

F r is the force exerted by the aluminium rod on the sample. It is determined by the gauges glued to the rod.

Gauges are also glued to the external face of the confinement tube so as to measure the axial profile of the circumferential strain (c.f. Figure 2 and Figure 3).

Sliding without fracture

A prior study of the competition between internal sliding [START_REF] Wiegand | Pressure and friction dependant mechanical strength -cracks and plastic flow[END_REF] and interfacial sliding is required to ensure the feasibility of this test. Indeed, further to a fracture, internal sliding of the sample may well be obtained in addition to the sliding at the interface. Internal sliding may even be obtained with some extrusion of material between the tube and the upper piston.

A sliding state (i) described by Coulomb's law between the I1 sample and the tube and a quasi-uniaxial strain state (ii) are considered. The sample is assumed to remain in the elastic domain. On the interface, the stress tensor is of the following form:

(3)

z r p fp p fp p , , 1 0 0 0 0 θ ν ν σ               - =
p being the interface pressure (positive compressive stresses), ν the Poisson coefficient of the material and f the friction coefficient at the interface. The expressions of the hydrostatic pressure P and the Von Mises stress σ mises are deduced as follows:

(4)

       +       - = + = 2 2 3 2 1 3 1 f p p P mises ν ν σ ν ν
The equation ( 2) must be satisfied if the sample is not to be fractured. The limit friction coefficient may thus be defined as:

(5) ( ) ( )           - -       + = 2 2 2 lim 2 1 3 1 3 1 ν ν α ν f
If f<f lim , the condition (2) is respected whatever the value of p. For f>f lim , the condition (2) is respected if p<p lim with p lim the limit pressure is defined by: (6)

ν ν α ν ν 3 1 3 2 1 2 2 lim + - +       - = f C p
p lim depends on ν and α, but a minor variation of ν may have a great influence.

Experimentally, the Poisson coefficient of the I1 is somewhere around 0.4. Figure 5 shows how this limit pressure varies as a function of ν and f. It is essential for the Poisson coefficient value ν to be known with suitable accuracy. It also shows that with a low Poisson coefficient value (0.38), pressures approaching 100MPa can be hoped for.

If f>0.45 is measured during the tests, it is almost certain that 100MPa will not be reached without fracturing the material. Naturally, in all cases, the samples must be checked a posteriori for any fractures.

Modelling of the test

The mechanical problem

The interface behaviour cannot be deduced by direct measurements. Indeed, a prior analysis of the test shows that the stress state at the interface is heterogeneous. A model must thus be established enabling this data to be acquired on the basis of measurements of F m and F r (c.f. Figure 6). For this, it is firstly necessary for the mechanical state of the sample to be determined.

Given that there is axial symmetry, the null components of the stress tensor σ(r,z) are those indicated in relation [START_REF] Azhdar | Determination of dynamic and sliding friction, and observation of stick-slip phenomenon on compacted polymers during high velocity compaction[END_REF]. The corresponding strains are also null. The displacement is axisymmetrical, the radial (oriented by u r ) and longitudinal (oriented by u z ) components depend on r and z.

The sample is forced into the tube by using the testing machine. Because of the difference in diameters u ri , the sample is subjected to pre-tightening linked to the radial displacement. If u r tube (r=R,z) represents the radial displacement withstood by the tube at r=R, then we obtain:

(7) ( ) ( ) ri tube r r u z R r u z R r u - = = = , ,
That is p(z) and τ(z) the normal and tangential stresses at the interface (at r=R): p(z)=σ rr (r=R,z) and τ(z)=σ rz (r=R,z). Conventionally, these are positive in compression.

Assuming that the friction at the interface is in accordance with Coulomb's law: the ratio of τ by p is equal to a coefficient f whose value is supposed to be little dependent on p.

However, the pressure range covered during the experiments is far-reaching and possible dependence of the friction coefficient f at pressure p must be envisaged [START_REF] Ben-Dor | Localized interaction models with non-constant friction for rigid penetrating impactors[END_REF]. The retained hypothesis is that a small variation of f as a function of p may be approached by an affine function:

(8) ( ) ( ) ( ) ( ) z p z p f z × = τ with ( ) ( ) p f p f β - = 1 0
Forces F m and F r are expressed in terms of boundary conditions (at z=0 and at z=L):

(9) ( ) ∫ = = R zz r dr z r r F 0 0 , 2 σ π (10) ( ) ∫ = = R zz m dr L z r r F 0 , 2 σ π
Given the tube's thickness and the gap between the rigidities of the materials in contact, the confinement tube is very rigid compared with the sample. In this case, no coupling between the behaviours of the two solids is considered and the presence of the tube may be modelled by a boundary condition of the imposed displacement type. According to [START_REF] Durand | Numerical solution of Cauchy problems in linear elasticity in axisymmetric situations[END_REF]:

(11) ( ) ( ) ri r tube r u z R r u z R r u - = = ⇒ = = , 0 ,
In the event of there being no excess thickness in the assembly and when the friction coefficient is not dependent upon the pressure, a two-dimensional approximate solution may be established.

(12) ( ) 0 , 0 = = ⇒ = z R r u u r ri
Assuming the sample is in an elastic state, the equilibrium equations, Hooke's law and the axial symmetry hypothesis are combined to produce the Navier equations:

(13) ( )( ) ( ) 0 2 1 1 2 , 2 , 2 , , 2 = - + + - + - zz r rz z r r r rr r u r u r u ru u r ν ν (14) ( ) ( )( ) 0 2 1 1 2 , , , , , = + - + + + - rr z r z rz r z r zz z ru u ru u ru ν ν
The friction law imposes at r=R:

(15) ( ) [ ] ( )( ) 0 2 1 1 2 , , , , = + - + + + - r z z r z z r r r u u r u r u u r f ν ν ν ν
Since the sample is confined, the radial displacement and its variations are very low with respect to the other terms and a first order approximation consists in leaving out certain terms to eventually obtain the following relations:

(16) ( ) 0 1 2 , , = + - rz z rr r u u ν (17) ( ) ( )( ) 0 2 1 1 2 , , , = + - + - rr z r z zz z ru u ru ν ν (18) ( ) 0 2 1 2 , , = - + r z z z u fu ν ν at R r =
A solution is sought after in the following form:

(19) ( ) ( ) ( ) ( ) ( ) ( )    + = = 0 , , u z r f z r u z r g z r u z r ϕ ψ
In the event of f being low, the above hypothesis is verified (c.f. appendix), and the following solution may be found:

(20) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )        +         - =       -         - - = ∫ ∫ 0 0 0 0 0 0 2 1 2 exp 0 1 2 exp 1 0 u r J R fz u d J d J R r R fz R f u z R r r γ ν ν ϕ ρ γρ ρ γρ ν ν ν ϕ ν
where J 0 is the zero order Bessel function of the first kind and γ is defined by:

(21) ( )( ) ν ν ν γ 2 1 1 2 2 - - = R f
The radial profile of the axial displacement u z (r) is, in fact, imposed by the relation [START_REF] Philippon | An experimental study of friction at high sliding velocities[END_REF]. The real boundary conditions at z=0 and at z=L may therefore not be respected exactly.

The form of this solution nevertheless enables the two-dimensional effects to be taken into account.

By combining the relations (20) with Hooke's law and with the boundary conditions ( 9)- [START_REF] Forquin | Influence of free water on the quasi-static and dynamic of strength of concrete in confined compression tests[END_REF], we obtain:

(22) ( )         - - = ν ν 1 2 exp R fL F F m r
This equation demonstrates that it is possible to evaluate f from the measurements of F m and F r when the geometry of the problem and the Poisson coefficient of the sample are known. The form of the stresses at the interface may also be obtained (c.f. appendix):

(23) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )                - - - =         - - - = ν ν π ν ν τ ν ν π ν ν 1 2 exp 1 1 2 exp 1 2 2 R L z f R fF z R L z f R F z p m m
This confirms the heterogeneity of the interface stresses mentioned above.

An approximate solution adapted to the test

The idea is now to simplify the model to obtain a solution that takes into account the effect of the initial excess of thickness of the sample and the dependency of the friction on the normal pressure.

The approach used is similar to that of Janssen [16]. It is based on two hypotheses: (i) the tube is assumed to be perfectly rigid and (ii) the axial, radial and circumferential stresses and strains do not depend on r. This enables the problem to be approached by a onedimensional incremental analysis, which is easier to perform.

The equilibrium equation ( 13) is no longer necessary and only the equilibrium equation ( 14) needs to be taken into account. We set: σ zz (r,z)=σ(z) and σ rr (r,z)=p(z). The equilibrium equation is thus:

(24) ( ) [ ] ( ) 0 , = + ∂ ∂ dz z d r r z r r rz σ σ
By integrating [START_REF] Vandersall | Experimental and modeling studies of crush, puncture, and perforation scenarios in the Steven impact test[END_REF] from r=0 to r=R, we obtain:

(25) ( ) ( ) dz z d R z σ τ 2 = and ( ) ( ) z R r z r rz τ σ - = ,
The components ε rr and ε θθ are positive in compression and assumed to be independent of r. Their values are:

(26) R u ri rr = = θθ ε ε
By applying Hooke's law, we then obtain:

(27) ( ) ( ) ( ) z z p R Eu ri νσ ν - - = 1
where E is the Young's modulus of the sample. It is equal to 2GPa.

The relations ( 9)- [START_REF] Forquin | Influence of free water on the quasi-static and dynamic of strength of concrete in confined compression tests[END_REF] become:

(28) ( )

0 2 = = z R F r σ π (29) ( ) L z R F m = = σ π 2
Because of the pre-tightening, a non-null force F m is required to make the assembly slide, even if F r is null. This force is written as F m0 . It is expressed as a function of the other parameters in the following form:

(30) ( ) ( ) ( )       -         - - +       - -       -         - = 1 1 2 exp 1 1 1 1 1 2 exp 0 0 0 ν ν ν β ν β ν π ν ν R L f R Eu Eu R Eu R L f F ri ri ri m
Let Q F and the constants Q F0 and r F F Q , be defined as:

(31) ( ) ( ) ( ) ( ) ( ) ( ) ( )                -                 -         - - +               - - - =                 -         - - +               - - = - = 2 0 0 , 2 0 0 0 0 1 1 1 2 exp 1 1 1 2 exp 1 1 1 2 exp 1 1 1 2 exp R R L f R Eu R L f Q R L f R Eu R L f Q F F F Q ri F F ri F m m r F r π ν ν ν ν β ν ν βν ν ν ν β ν ν
By combining the equilibrium equation, the laws of friction, Hooke's law and the boundary conditions, we obtain the following affine relation between Q F and F r :

(32) r F F F F F Q Q Q r × + = , 0
By adding the boundary conditions, we finally obtain [START_REF] Rajagopalan | Novel experimental techniques for investigating time resolved high speed friction[END_REF]. This shows that the radial independence hypothesis with respect to the axial, radial and circumferential stresses and strains has no incidence on the simultaneous evolutions of F m and F r and thus no incidence on the processing method.

The evolutions of F m and F r obtained during the tests have enabled three unknowns f 0 , β and u ri to be identified. Only f 0 and β are really of interest since these are the parameters which control the behaviour of the interface. However, it is important for u ri to be known thereafter in order to be able to accurately simulate the test.

Lastly, the normal and tangential stresses at the interface are determined using the following formulae:

(33) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )                        + - +       -         - -         - -       - + - = + - +       -         - -         - - = 2 2 0 0 2 0 2 0 0 1 1 1 2 exp 1 2 exp 1 1 1 1 2 exp 1 2 exp m ri m ri m ri F REu R R L z f R L z f F REu R f z F REu R R L z f R L z f z p ν π π ν ν ν β ν ν β ν π π ν τ ν π π ν ν ν β ν ν
Similarly, when the boundary conditions are introduced into the expression of the interface stresses, equality (23) is obtained.

Experimental results

Performed tests

Two types of tests were performed, during one of which the strains along the tube were measured by the eight circumferential gauges (c.f. Figure 2 and Figure 3). For practical reasons, this test was performed at low pressure. The other test was performed up to a high pressure of 100MPa. Hereafter, the first test will be termed "lower pressure test" and the second "higher pressure test". The sample is merely pushed on the spring by the testing machine (c.f. Figure 2) as the measurements are taken. Thus, the normal pressure (by Poisson effect on the sample) and the sliding are imposed simultaneously thanks to the spring stiffness (c.f. section 2.1). Both tests were carried out with the same tube and the same sample.

The measurements allow the values of the forces F m and F r to be accessed during loading. The raw results for both tests are presented in Figure 7. I1 properties do not vary a lot from a sample to another and friction tests are therefore reproducible. These data are then typical of several measurements.

The idea is to analytically link (using relations (30)-( 31)-( 32)) the friction parameters and other unknowns to magnitudes identifiable from the evolution of the two forces. This, firstly, enables the coherence of the analytical model to be verified.

An order of f relative uncertainty

Analysis of the lower pressure test

Given that F m0 =2.81kN is obtained by the linear extrapolation of the blue curve in Figure 7, the evolution of Q F is drawn as a function of F r (c.f. Figure 8).

The Q F ratio seems constant and independent of F r . It can thus be assimilated to its mean value: Q F =0.41. Eliminating the dependency of this ratio to F r comes down to equating β=0 in (32), which becomes:

(35) ( ) ( ) F Q L R f ln 2 1 0 ν ν - - =
(30) can thus be written as:

(36) ( ) F F m ri Q RE Q F u - = 1 0 π ν
There are still two parameters remaining to be identified: u ri and f 0 . Since the other parameters are known, identification is made using the measured values of Q F and F m0 . We deduce from ( 35) and (36) that u ri =0.0124mm and f 0 =0.167.

Analysis of the higher pressure test

For this test, the first measuring point corresponds to F m0 and gives F m0 =900N. The experimental evolution of Q F may thus be drawn as a function of F r (c.f. Figure 9).

The relations (30)-( 31)-( 32) enable the identification of the parameters f 0 , β and u ri based on the experimental values of the slope r F F Q , and the intercept Q F0 of the linear regression line and on the experimental value of F m0 . Since the system is difficult to be solved, an approximation is added: the initial force F m0 being relatively low with respect to the values attained by the forces F m and F r , u ri is neglected with respect to the other terms in (32). The following equations are thus obtained:

(37) ( ) ( ) 0 0 ln 2 1 F Q L R f ν ν - - = (38) ( ) ( ) 0 , 2 1 1 F r F F Q Q R - - = ν π ν β
A linear regression leads to Q F0 =0.339 and r F F Q , =0.00223kN -1 . Since the other parameters are known, f 0 =0.203 and β=0.00159MPa -1 can be deduced.

In relation (30), the second order terms at u ri are neglected and by taking relation (37) into account, we obtain:

(39) ( ) 0 0 0 1 F F m ri Q RE Q F u - = π ν
(39) gives u ri =0.00294mm.

The system of relations (33) and the numerical simulations (c.f. section 5.1) enable the contact pressures to be estimated at the interface, whereas these are not able to be measured. It is thus possible to account for the pressure ranges covered during this test (c.f. Figure 10).

The friction coefficients identified during the two tests are not quite equal. At a pressure of 30MPa, the first test gives a coefficient of 0.17 and the second one a coefficient of 0.19. The interface conditions change over time (polishing of the sample, smoothing of any asperities of the steel, etc) since these tests are performed successively with the same tube and the same sample, the assembly being simply put back into position after each test. The lower pressure test being performed after the higher pressure test, the reduction in friction between the two tests can be explained by the changing interface conditions. Moreover, the sliding rates imposed during the lower pressure test are of a magnitude of a mm/h whereas they are of a magnitude of a mm/min for the high pressure test. This change may also have an effect.

The value of the friction coefficient reduces when the contact pressure increases, as suggested by certain authors [START_REF] Ben-Dor | Localized interaction models with non-constant friction for rigid penetrating impactors[END_REF]. For the example shown here, this value goes from 0.2 at very low contact pressure to 0.16 at contact pressure of 140MPa.

Simulations and discussion

The tests were numerically simulated to validate the analysis and the friction parameters thus identified.

Modelling was performed using the finite element method. The software used was ABAQUS CAE / Standard (implicit). The performed calculations were two-dimensional and axisymmetric. Quadrangular elements with quadratic interpolation were used. All the normal contacts were defined by "hard contact" "direct contact" (no interpenetration), separation being allowed (no bonding). The tangential contact between the confinement tube and the sample obeys a Coulomb's law (the coefficient being an affine function of the pressure) imposed using the method of Lagrange multipliers. The other tangential contacts are without friction.

The mechanical properties selected for the steel and the aluminium were respectively a Young's modulus of 210GPa and a Poisson coefficient of 0.33, and a Young's modulus of 74GPa and a Poisson coefficient of 0.3. The sample behaviour is assumed to be perfectly elastic with a 2GPa Young's modulus (c.f. section 3.2) and a 0.4 Poisson coefficient (c.f. section 2.2). This assumption will be checked in section 5.4.

The finite elements are squares of 1mm in length. There are therefore 10 elements in the radius of the sample.

Validity of the analysis method

The test was simulated with the parameters determined during the analysis. Numerical simulation enables the evolutions of F m and F r (and thus of Q F ) to be obtained as well as the interface stress profiles (which is to say p(z) and τ(z)). Thereafter, the results were compared with those obtained experimentally.

For the lower pressure test, loading was simulated with F m rising up to F m MAX =17kN

(which is equal to the maximal F m experimentally reached). The parameters identified from the model allow the finite element evolution of Q F as a function of F r to be very close to the experimental one (c.f. Figure 8). The friction parameter identification technique is thus validated. Figure 11 shows that our approach cannot be used to predict boundary effects, but this has no incidence on the identification of the friction parameters.

For the higher pressure test, loading was simulated with F m rising up to F m MAX =75kN

(which is equal to the maximal F m experimentally reached). The conclusions regarding the numerical simulation of the higher pressure test are exactly the same as previously (c.f. Figure 9 and Figure 11).

Figure 12 displays (for the higher pressure test) that the radial pressure field is homogenous in the radial direction everywhere in the sample but near the top and bottom boundaries.

Despite its imperfections, the model enables the simultaneous evolutions of the magnitudes Q F and F r to be accurately obtained, and this for both tests. Q F and F m0 are both calculated with f 0 =0.167, β=0 and u ri =0.0124mm using relations (30)-( 31)-(32), using the simulation of the lower pressure test and using an additional simulation considering a rigid tube. The results (c.f. Table 1) show that the slight deviations between the model and the simulations are actually due to the flexibility of the tube which is not taken into account in the model. Indeed, the relative deviations between the model and the simulation are higher than those between the model and the simulation considering a rigid tube.

The numerical simulations enable the validity of the analysis, and thus the validity of the friction parameters thereby defined, to be ensured.

Consistency of superabundant measurements

During a test, the circumferential strain profile along the tube ε θθ (z) was measured at several steps during the loading, that is to say at several values of the force F m . A discrete experimental profile was obtained and compared with the numerically obtained profile so as to verify their consistency. Figure 13 shows the position and width of each of the gauges as well as their measurements.

With the exception of the values obtained at z=L, a satisfactory correspondence is to be noted between the experimental and numerical axial profiles (c.f. Figure 13).

The higher pressure test was performed at a higher loading rate. For practical reasons, only the strain of the gauge glued at z=0.5L was measured (c.f. Figure 14).

Despite some imperfections, no significant deviations were observed between the simulated strains and the measured ones, thereby providing an additional guarantee with regard to the consistency of the method.

Influence of the friction at the top and bottom boundaries

The interface conditions at the top and bottom boundaries (located at z=0 and at z=L) may have an influence on the behaviour of the sample, in particular during the tests on the Hopkinson bars [START_REF] Hartley | An experimental assessment of friction effects in the split Hopkinson pressure bar using the ring compression test[END_REF], [START_REF] Li | Further investigation on the dynamic compressive strength enhancement of concrete-like materials based on split Hopkinson pressure bar tests[END_REF], [START_REF] Walley | A high-speed photographic study of the rapid deformation of metal annuli: comparison of theory with experiment[END_REF]. It is thus necessary to ensure that these effects remain negligible with the selected configuration (slender sample with radial confinement).

For this, the numerical results obtained for two extreme cases were compared. The results presented in sections 5.1 and 5.2 were obtained by modelling the contacts at the top and bottom boundaries without friction (the friction coefficient being assumed to be null). The simulations were performed again but with changed boundary conditions: instead of imposing a null friction stress at the boundaries, contacts with null relative displacement were imposed (friction coefficient being assumed to be infinite).

We observed that the friction coefficient value at the top and bottom boundaries has no influence (relative deviation of less than 2%) on the evolution of Q F as a function of F r , and thus no influence on the parameters identified during the analysis.

Verification of the elasticity hypothesis

The elasticity hypothesis of the sample had to be verified to ensure the consistency of the analysis. The criterion (σ mises -αP) was calculated in the sample during the numerical simulations and this after the loading when the forces were maximal. Thereafter, the fulfilment of the condition (σ mises -αP)<C (where α=0.64 and C=25MPa) is verified.

The numerical results obtained for the two extreme cases at the boundaries (null friction coefficient, infinite friction coefficient) were compared. These calculations were performed only for the higher pressure test for which the stresses are the highest.

The value of the criterion was at most 13MPa; this remains therefore within the elastic domain (c.f. Figure 16).

The simulations were performed once again but by imposing an infinite friction coefficient at the top and bottom boundaries. Loading was also simulated with F m rising up to

F m MAX =75kN.
The value of the criterion was at most 5.1MPa; this remains therefore within the elastic domain (c.f. Figure 17).

The values attained by the criterion (σ mises -αP) depend on the contact conditions at the boundaries. In the two envisaged extreme cases, we observed that the elastic limit C was never attained. Our analysis, which uses an elastic solution, is therefore valid. Significantly, the criterion (σ mises -αP) is only positive in very localized areas.

The sample is forced into the tube by using the testing machine. Numerical simulations of this insertion display that (2) remains satisfied for the lower and the higher pressure test. As a result, no plastic deformations take place during insertion.

This article presents a new experimental test for friction between steel and a brittle material. A technique from the compacting of granular materials is used so as to maintain the integrity of the brittle material under high interface pressures.

The analysis, performed analytically and numerically, enables the friction parameters to be obtained and shows that pressures of around 100MPa are reached. The simulations also show that the sample does not fracture under the load to which it is subjected. The developed assembly thus enables pressures with higher magnitudes than in [START_REF] Dickson | Frictional Heating and Ignition of Energetic Materials[END_REF] and [START_REF] Hoffman | Aspect of the tribology of the plastic bonded explosive LX-04[END_REF] to be attained without any deterioration of the material. The sliding rates imposed at the interface are of an order of magnitude of a mm/min (higher pressure test) and a mm/h (lower pressure test). The measured friction coefficient is of around 0.2. The device enables the variation of this coefficient to be quantified when the contact pressure increases.

The device enables something to be made to slide under high contact pressure, and the method has been validated. Test campaigns may be undertaken on inert or reactive materials.

The dispersion aspect of the parameters identified may be understood.

One alternative to the analysis technique would be the used of an inverse method. The feasibility of such a method is shown in [START_REF] Durand | Numerical solution of Cauchy problems in linear elasticity in axisymmetric situations[END_REF].

The next stage of the study is to mount the experimental device on a system of Hopkinson bars so as to reach sliding rates of an order of magnitude of 10m/s. Result of the finite element simulation (red) detailed in section 5. Result of the numerical simulation detailed in section 5 (red).

∆FF

  by introducing u ri =0 and β=0 in relations (30)-(31)-(32), which leads to: ratio is around 0.4 for both tests. F m is measured by the testing machine sensors (accuracy approx. 1N) and F r is deduced from strain gauges (c.f. section 2.1). As the machine sensors are very accurate

Figure 1 :

 1 Figure 1: Numerical simulation of the heating of an explosive projected at 100m/s against a

Figure 2 :

 2 Figure 2: Diagram of the device placed in a testing machine. The sample is 20mm in diameter

Figure 3 :

 3 Figure 3: Photograph of the device placed in a testing machine.

Figure 4 :

 4 Figure 4: Definition of the forces on the assembly.

Figure 5 :

 5 Figure 5: Evolution of p lim as a function of f for different values of ν.

Figure 6 :

 6 Figure 6: Definition of the mechanical problem.

Figure 7 :

 7 Figure 7: Comparison of the evolutions of F r as a function of F m during the two tests.

Figure 8 :

 8 Figure 8: Experimental evolution of Q F as a function of F r for the lower pressure test (blue).

Figure 9 :

 9 Figure 9: Experimental evolution of Q F as a function of F r for the higher pressure test (blue).

Figure 10 :

 10 Figure 10: Friction law identified using the higher pressure test (pressure range between 5

Figure 11 :

 11 Figure 11: Comparison of the analytical and numerical profiles of pressure p(z) for both tests when the force is maximal (F m =F m MAX ).

Figure 12 :

 12 Figure 12: Value of the radial pressure p in MPa for the higher pressure test.

Figure 13 :

 13 Figure 13: Comparison of the experimental and numerical profiles of ε θθ (z) for two values of

Figure 14 :

 14 Figure 14: Comparison of the experimental and numerical evolutions of ε θθ (z=0.5L) as a

Figure 15 :

 15 Figure 15: Comparison of the experimental and numerical evolutions of ε θθ (z=0.5L) as a

Figure 16 :

 16 Figure 16: Value of the criterion (σ mises -αP) in MPa after the loading for the higher pressure

Figure 17 :

 17 Figure 17: Value of criterion (σ mises -αP) in MPa after the loading for the higher pressure test

Figure 18 :

 18 Figure 18: Evolution of the relative error err as a function of the friction coefficient f.

Table 1 :

 1 Q F0 and F m0 calculated with f 0 =0.167, β=0 and u ri =0.0124mm and relative deviations between analytical relations and simulations.

	Q F (relative deviation) 0.410	0.408 (-0.49%)	0.421 (2.7%)
	F m0 (relative deviation) 2.80kN	2.83kN (1.1%)	2.68kN (-4.3%)

Appendix

By using the form of u z defined in relation [START_REF] Philippon | An experimental study of friction at high sliding velocities[END_REF], relation [START_REF] Li | Further investigation on the dynamic compressive strength enhancement of concrete-like materials based on split Hopkinson pressure bar tests[END_REF] produces the solution:

and then, relation [START_REF] Kim | Nanostructures generated by explosively driven friction: Experiments and molecular dynamics simulations[END_REF] becomes:

The form of f(r) is deduced from (A2) and this gives:

where J 0 is the zero order Bessel function of the first kind.

The relation (16) enables u r (r,z) to be determined in the form defined by [START_REF] Philippon | An experimental study of friction at high sliding velocities[END_REF]. Taking [START_REF] Hartley | An experimental assessment of friction effects in the split Hopkinson pressure bar using the ring compression test[END_REF] into account, we then obtain:

The equation (A3) enables the value of γ to be determined. To solve this equation, a hypothesis is first posed, then verified a posteriori: the terms at (γR) n where n≥3 are assumed to be negligible and the Bessel functions at γR and at γr can be assimilated to their limited development of the second order. In this case we obtain:

where J 1 is the first order Bessel function of the first kind.

The equation (A3) is thus easily solved, and we obtain:

To ensure the consistency of the approach, the accuracy of the approximation (A7) must be evaluated. For this, we must calculate the relative error err committed as a function of the value of the friction coefficient f.

By substituting the expression given by (A8), we obtain:

According to Figure 18, the approximation (A7) remains valid for the usual values of the friction coefficient f. The formulae [START_REF] Rajagopalan | An experimental method to study high speed sliding characteristics during forward and reverse slip[END_REF] are determined on the basis of Hooke's law, neglecting u r and u r,r and applying the approximation (A7).