Bernhard Beckermann

George Labahn
email: glabahn@uwaterloo.ca

Ana C Matos
email: matos@math.univ-lille1.fr.

On rational functions without Froissart doublets

Keywords: numerical analysis, rational functions, Padé approximation, Froissart doublets, spurious poles, numerical coprimeness AMS subject classification: 41A21, 65F22

In this paper we consider the problem of working with rational functions in a numeric environment. A particular problem when modeling with such functions is the existence of Froissart doublets, where a zero is close to a pole. We discuss three different parameters which allow one to monitor the absence of Froissart doublets for a given general rational function. These include the euclidean condition number of an underlying Sylvester-type matrix, a parameter for determining coprimeness of two numerical polynomials and bounds on the spherical derivative. We show that our parameters sharpen those found in a previous paper by two of the authors.

Introduction

Let C[z] be the space of polynomials with complex coefficients, C n [z] the subset of polynomials of degree at most n, and

C m,n [z] = p q : p ∈ C m [z], q ∈ C n [z], q = 0
the set of rational functions. Rational functions have long played an important role in applied mathematics. As an example, Padé approximants and rational interpolants are used for approximation, analytic continuation and for determining singularities of a function [START_REF] Baker | Padé Approximants[END_REF][START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]. Also, Padé approximants of the z transform of noisy signals are employed for detecting the number of significant signals, their frequencies, damping, phase and amplitude [START_REF] Bessis | Padé approximations in noise filtering[END_REF]. Other applications include sparse interpolation [START_REF] Giesbrecht | Symbolic-numeric Sparse Interpolation of Multivariate Polynomials[END_REF][START_REF] Kaltofen | On exact and approximate interpolation of sparse rational functions[END_REF], computer algebra [START_REF] Zur Gathen | Modern Computer Algebra[END_REF][START_REF] Driver | Padé approximations to the Logarithm II: Identities, Recurrence and Symbolic Computation[END_REF] and exponential analysis [START_REF] Cuyt | Sparse interpolation and Rational approximation[END_REF][START_REF] Peelman | Sparse reconstruction of correlated multichannel activity[END_REF][START_REF] Potts | Parameter estimation for nonincreasing exponential sums by Prony-like methods[END_REF].

In order to successfully use rational functions for modelling, one first has to address the subtle question of choosing a priori the degrees m, n. We want to make sure that 1 the rational function r ∈ C m,n [z] is nondegenerate, that is, at least one of the numerator or denominator degrees is maximal after removing common factors. In addition, such a rational function should also be sufficiently "far" from C m-1,n-1 [z]. Indeed, overshooting the degree may produce strange artifacts commonly referred to as spurious poles. By this we mean we might have so-called Froissart doublets [START_REF] Froissart | Approximation de Padé: application à la physique des particules élémentaires[END_REF], that is, a pair of points, one a pole and the other a zero of r, which are close to each other. This of course makes it impossible to approach smooth functions with such rational functions. Another second well known artifact is a simple pole with small residual, which does not seem to be significant if one wants to evaluate r at points not too close to such a pole. These issues become particularly significant when computation is done in a numeric environment where one can only obtain close rather than exact answers.

In a recent paper [START_REF] Gonnet | Robust Padé approximation via SVD[END_REF], the authors introduced the notion of robust Padé approximants, a lower order Padé approximant based on the SVD of the underlying Toeplitz matrix. There the authors showed in many illustrating examples that their robust Padé approximants no longer have spurious poles. Only later was it shown that the underlying nonlinear Padé map taking the coefficients of the initial Taylor series and mapping them to the coefficients in the basis of monomials of the numerator and denominator of a Padé approximant is forward well conditioned (but not necessarily backward) at such robust Padé approximants [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Theorem 1.2], and that robust Padé approximants may have spurious poles [START_REF] Mascarenhas | Robust Padé Approximants may have spurious poles[END_REF][START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]. A first important contribution to this set of questions is [START_REF] Werner | On the continuity of the Padé operator[END_REF] for the continuity of the Padé map.

In computer algebra the area of symbolic/numeric computation often considers correctness and stability issues when working with polynomial arithmetic. For example there has been considerable work on problems such as the numerical gcd of two polynomials having floating point coefficients. However there seems to be very little work which deals with numerical analysis around rational functions. Indeed even the first issue of clarifying how to measure distances in C m,n [z] has not really been considered.

It seems natural that one should expect a connection between "nearly" degenerate rational functions and numerators and denominators having a non-trivial numerical gcd. This includes the two papers [START_REF] Beckermann | When are two numerical polynomials relatively prime[END_REF][START_REF] Corless | The singular value decomposition for polynomial systems[END_REF] where coprimeness parameters are considered and which both make the link with the underlying Sylvester matrix formed by the coefficients of the numerator and denominator (see Definition 2.1).

The aim of this paper is to discuss three different parameters which allow one to monitor the absence of Froissart doublets for a given general rational function r ∈ C m,n [z].

• The euclidean condition number of underlying Sylvester-type matrices depending on some integer ;

• the coprimeness parameter of [START_REF] Beckermann | When are two numerical polynomials relatively prime[END_REF][START_REF] Corless | The singular value decomposition for polynomial systems[END_REF];

• bounds on the spherical derivative.

In each case we will show how our first two parameters generalize and sharpen the parameters presented in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]. In the case of the spherical derivative, the two new parameters introduced here are essentially best Lipschitz constants of a rational function and we show how measuring distances in C m,n [z] also partially sharpens some distance measures found in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF].

This paper is a follow up to the paper [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF] where some of the same problems were considered. In order to explain our contributions in more detail, we describe results from the previous paper along with our new findings. For this it is helpful to refer to different properties given in Figure 1. The authors in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF] refer to r = p/q with modest euclidean condition number of the underlying Sylvester-type matrix for = 1 as well conditioned rational functions, and deduced several properties of such functions. This includes, for example, the absence of Froissart doublets [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Theorem 1.3(a)] and of small residuals for simple poles [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Theorem 1.3(b)], a large distance to the set C m-1,n-1 [z] of degenerate rational functions [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Theorem 1.4], but also a modest forward and backward condition number for the non-linear Padé map for well conditioned Padé approximants [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Theorem 1.2]. They also establish the equivalence of two different distances in C m,n [z], one based on values and the other on coefficients of rational functions [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Theorem 4.1]. This corresponds to the implications on the right of Figure 1. We will establish later in Theorem 2.3 that the choice of our parameter in Definition 2.1 is not essential.

The results in this paper correspond to the left side of Figure 1. In §2.1 we recall some of the findings from [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF] and also the coprimeness parameter of [START_REF] Beckermann | When are two numerical polynomials relatively prime[END_REF][START_REF] Corless | The singular value decomposition for polynomial systems[END_REF]. We show, in Theorem 2.5, how this coprimeness parameter allows one to monitor the absence of Froissart doublets, and so generalizes and sharpens the previous attempts found in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF][START_REF] Beckermann | When are two numerical polynomials relatively prime[END_REF]. In §2.2 we introduce two new parameters based on the spherical derivative. We show, in Theorem 2.9(a),(b) and Corollary 2.8, that these parameters are essentially best Lipschitz constants of a rational function, and that these new parameters also allow one to insure the absence of Froissart doublets and poles with small residual. In addition, in Theorem 2.12 we describe some special cases where these findings are sharper than those of Theorem 2.5. Finally, in §2.3 we come back to the question of comparing distances in C m,n [z]. We show, in Theorem 2.13, that [4, Theorem 4.1] can be partly sharpened in terms of the coprimeness parameter and give an example showing that a second inequality cannot be improved. This completes the picture of Figure 1.

The remainder of this paper is as follows. Statements of our main results involving our three parameters are presented in three subsections in §2, with proofs of all our main statements given in §3. The paper then ends with a conclusion and topics for future research in §4.

Main results

In this section we present the main results mentioned in the previous section. Here all theorems are stated with the proofs given later in the following section.

Measure of coprimeness and Sylvester type matrices

In what follows we consider fixed integers m, n ≥ 0. In order to simplify notation, we will not explicitly indicate the dependency on m, n of each object. For a polynomial c(z) = c 0 + c 1 z + • • • + c n z n with coefficients c j we denote by vec(c) = (c 0 , c 1 , • • • , c n ) T its coefficient vector, with the size of this vector being clear from the context. We start by introducing a so-called Sylvester-type matrix S ( ) associated to a pair of polynomials.

∈ C m [z], q ∈ C n [z]
, with coefficients p j , q j , respectively, the associated (m + n + ) × (m + n + 2 ) Sylvester-type matrix of p and q is defined by

S ( ) (p, q) =               p 0 p 1 p 0 . . . p 1 . . . p m . . . . . . p 0 p m p 1 . . . . . . p m n+ q 0 q 1 q 0 . . . q 1 . . . q n . . . . . . p 0 q n q 1 . . . . . . q n m+               ∈ C (m+n+ )×(m+n+2 ) .
When = 0, S ( ) (p, q) reduces to the transpose of the classical Sylvester matrix [START_REF] Geddes | Algorithms for Computer Algebra[END_REF] while when = 1 we get the Sylvester-type matrix used in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]. The more general allows us to consider increased degrees in an associated diophantine equation connected to polynomial gcd computation of p and q.

It is well known [START_REF] Geddes | Algorithms for Computer Algebra[END_REF] that the classical square Sylvester matrix S (0) (p, q) is invertible if and only if the polynomials p and q are coprime and the defect min(m-deg (p), n-deg (q)) is equal to zero, that is, the rational function p/q is nondegenerate. More generally, S ( ) (p, q) has full row rank if and only p/q is nondegenerate. We refer to [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Lemma 3.1] for a proof in the case = 1, while a proof for > 1 is similar, based on the relation 1, z, ..., z m+n+ -1 S ( ) (p, q)

(2.1)

= p(z), z 1 p(z), ..., z n+ -1 p(z), q(z), z 1 q(z), ..., z m+ -1 q(z) .

In order to make the link with the coprimeness parameter discussed by Corless, Gianni, Trager, and Watt in [START_REF] Corless | The singular value decomposition for polynomial systems[END_REF] we introduce as in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]

a norm in C[z] × C[z] through the formula (p, q) 2 = vec(p) 2 2 + vec(q) 2 2 ,
and consider the following quantities.

Definition 2.2. For p ∈ C m [z], q ∈ C n [z], and a set K ⊆ C, consider 2 (p, q) = inf (p -p, q -q) 2 : ( p, q) ∈ C m [z] × C n [z] have a common root , K,2 (p, q) = inf z∈K |p(z)| 2 m j=0 |z| 2j + |q(z)| 2 n j=0 |z| 2j 1/2 .
The coprimeness parameter 2 (p, q) measures the distance to the set of pairs of polynomials with a non-trivial gcd. At the same time it is mentioned in [6, Remark 4] that 2 (p, q) coincides with C,2 (p, q), with the latter quantity being much more accessible since one minimizes only with respect to the single complex parameter z.

Since S ( ) (p, q) 2 is not too far from (p, q) 2 , our coprime measure 2 (p, q) approximately also gives the distance of S ( ) (p, q) to the set of singular Sylvester-type matrices S ( ) ( p, q), that is, a kind of smallest structured singular value [START_REF] Rump | Structured Perturbations Part 1: Normwise Distances[END_REF]. As such, a bound of the form

C,2 (p, q) = 2 (p, q) ≥ 1 √ m + n + 1 S ( ) (p, q) † 2 (2.2)
in terms of the norm of the pseudo-inverse is not surprising. To see this we just take norms in (2.1) (see also [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Lemma 5.1] for a proof in the case = 1). In [4, §6.2] we conjectured that the dependency on of the right-hand side of (2.2) is not important. In the present paper we are able to state:

Theorem 2.3. Let p ∈ C m [z], q ∈ C n [z] be such that p/q is nondegenerate. Then S (0) (p, q) -1 2 ≤ S ( ) (p, q) † 2 ≤ (1 + √ ) S (0) (p, q) -1 2 (2.3)
for all integers ≥ 0.

Remark 2.4. The authors in [START_REF] Beckermann | When are two numerical polynomials relatively prime[END_REF] have obtained more compact expressions by choosing in Definition 2.2 a different norm for pairs of polynomials, namely (p, q) 1 = max( vec(p) 1 , vec(q) 1 ).

This allowed them to deduce that S ( ) (p, q) 1 = (p, q) 1 , independent of . In this case, the one-norm equivalent of Definition 2.2 becomes It was also shown in [START_REF] Beckermann | When are two numerical polynomials relatively prime[END_REF]Theorem 4.1] that 1 (p, q) = C,1 (p, q), where

K,1 (p, q) := inf z∈K max |p(z)| max (1, |z| m ) , |q(z)| max (1, |z| n ) = inf z∈K (1, z, • • • , z m+n+ -1 )S ( ) (p, q) 1 (1, z, • • • , z m+n+ -1 ) 1 .
The last relation implies 1 (p, q) ≥ 1 S (0) (p,q) -1 1 , as mentioned already in [3, Lemma 2.1].

Froissart doublets

Let us now turn to the question of existence of a Froissart doublet for a rational function

r = p/q ∈ C m,n [z],
that is, a pair consisting of a zero z p and a pole z q of r which are close to each other. In [4, Theorem 1.3(a)] it was shown that

|z p -z q | ≥ 1 3 √ 2(m + n + 1) 3/2 cond(S (1) (p, q)) , (2.4)
provided that both z p and z q are in the closed unit disk D. Here cond(B) = B 2 B † 2 denotes the condition number with respect to the euclidean norm. It seems reasonable to expect that a sufficiently large s (p, q) for s ∈ {1, 2} also implies the absence of Froissart doublets, since in this case p, q is relatively far from a pair of polynomials having a non-trivial gcd. Let χ(x, y) := |x -y|

1 + |x| 2 1 + |y| 2 .
be the chordal metric obtained by taking the euclidean distance on the Riemann sphere S 2 which is identified with the extended complex plane C ∪ {∞} through stereographic projection. Then the distance between a pole and a zero of a rational function as measured using the chordal metric is bounded from below by:

Theorem 2.5. Let K ⊆ C and r = p q ∈ C m,n [z].
Then for any pair z p , z q ∈ C with p(z p ) = 0, q(z q ) = 0 and s ∈ {1, 2} we have that

K,s (p, q) 2 max ( m vec(p) s , n vec(q) s ) ≤ χ(z p , z q ). (2.5)
Moreover, if K ⊆ D or 1/K ⊆ D we can replace the maximum in the denominator by a minimum.

It follows from Theorem 2.5 that if p and q are numerically coprime (that is, having a sufficiently large K,s (p, q)) then r = p/q cannot have any Froissart doublets. Notice that (2.5) for s = 2 and K = C is sharper than (2.4). Indeed using (2.5) one can obtain

|z p -z q | ≥ 1 2(m + n + 1) 3/2 cond(S (1) (p, q)) ,
by verifying the two inequalities max(m vec(p) 2 , n vec(q) 2 ) ≤ (m + n + 1) S (1) (p, q) 2 and |z p -z q | ≥ χ(z p , z q ), and by applying the bound (2.2). Special cases of Theorem 2.5 have been claimed without proof in [3, §4] for s = 1 and K = C, and established in [4, Lemma 6.1] for s = 2 and K = D.

It is interesting to note that the indicators used in (2.4) and (2.5) are not sensitive with respect to a small perturbation of the numerator and denominator. Here we can state the following.

Theorem 2.6. Let K ⊆ C and p q , p q ∈ C m,n [z]. (a) If p
q is nondegenerate and

(p -p, q -q) 2 ≤ 1 3 √ m + n + 1 S (1) (p, q) † 2 then 1 2 ≤ cond(S (1) ( p, q))/cond(S (1) (p, q)) ≤ 2. (b) Let s ∈ {1, 2}. If (p -p, q -q) s ≤ 1 2 K,s (p, q) then 1 2 ≤ K,s ( p, q)/ K,s (p, q) ≤ 3 2 .
Notice that, according to (2.2), the neighborhood in part (b) for s = 2 is larger than the neighborhood in part (a).

We are now able to show that inequalities (2.4) and (2.5) are robust in the sense that they remain valid up to some modest constant if z p and z q are roots not of (p, q) but of some perturbed ( p, q) ∈ C m [z] × C n [z] sufficiently close to (p, q). For (2.4) this has been done before in [4, Theorem 1.3(a)], and we essentially repeat their arguments. For (2.5), it is convenient to write first the slightly weaker inequality

χ(z p , z q ) ≥ K,s (p, q) 2(m + n) (p, q) s ,
and to observe that K,s (p, q) ≤ (p, q) s . Then Theorem 2.6 yields the following.

Corollary 2.7. Let K ⊆ C and p q , p q ∈ C m,n [z] and z p , z q ∈ K with p(z p ) = q(z q ) = 0. Then, under the assumptions of Theorem 2.6(a), 1) (p, q)) .

|z p -z q | ≥ 1 6 √ 2(m + n + 1) 3/2 cond(S (
Furthermore, under the assumptions of Theorem 2.6(b),

χ(z p , z q ) ≥ K,s (p, q) 6(m + n) (p, q) s .
Again, using (2.2) one may show that the second statement for s = 2 implies the first one.

Spherical derivatives and small residuals

In this subsection we will introduce a new parameter in order to monitor the existence of Froissart doublets. Recall that the spherical derivative of a rational function r ∈ C m,n [z] is given by

ρ(r)(z) = |r (z)| 1 + |r(z)| 2 (2.6)
while for any K ⊆ C we set ρ K (r) := sup z∈K ρ(r)(z).

(2.7)

Note that 1 ρ(r)(zq) equals the modulus of the residual of a simple pole z q . Hence the following statement, complementing [4, Theorem 1.3(b)], is immediate.

Corollary 2.8. Let β be the residual of a simple pole z q of r in K. Then |β| ≥1 ρ K (r) .

We show below that the quantity ρ K (r) is the best Lipschitz constant for r in K. As such, for a reliable evaluation of r(z) for z ∈ K, it seems to be reasonable to restrict ourselves to rational functions r with modest ρ K (r). On the other hand, if we want to measure the distance of arguments in terms of the chordal metric, another indicator is more appropriate, namely

ν(r)(z) = (1 + |z| 2 ) |r (z)| 1 + |r(z)| 2 and ν K (r) := sup z∈K ν(r)(z). (2.8)
Let us now turn to Froissart doublets and compare our new indicators with those given previously.

Theorem 2.9. Let K ⊆ C and r = p q ∈ C m,n [z] with p and q coprime and z p , z q ∈ K with p(z p ) = q(z q ) = 0.

(a) If K is convex then ρ K (r) = sup z 1 ,z 2 ∈K χ(r(z 1 ), r(z 2 )) |z 1 -z 2 | . (2.9)
In particular,

|z p -z q | ≥ 1 ρ K (r) . (b) If K is spherically convex 1 then ν K (r) ≤ sup z 1 ,z 2 ∈K χ(r(z 1 ), r(z 2 )) χ(z 1 , z 2 ) ≤ π 2 ν K (r). ( 2

.10)

In particular, χ(z p , z q ) ≥2 πν K (r) . It is also interesting to explore the links between the spherical derivative and the numerical measure of coprimeness. Here the following observation is helpful.

Remark 2.10. Notice that, by definition K,s (p m , q m ) = K,s (p, q) m strongly depends on m, whereas ρ K (r m )(z) ≤ 2mρ K (r)(z) for any rational function r. This follows from

ρ(r m )(z) = m |r(z)| m-1 |r (z)| 1 + |r(z)| 2m ≤ |r(z)| m-1 max(1, |r(z)| 2m-2 ) m |r (z)| max(1, |r(z)| 2 ) ≤ m 2 |r (z)| 1 + |r(z)| 2 = 2mρ(r)(z),
using the fact that 2 max(1, |r(z)| 2 ) ≥ 1 + |r(z)| 2 and then taking sup over K.

In the following example, which was already studied in [3, Example 5.3], we see that the bound of Theorem 2.9(a) is approximately sharp whereas Theorem 2.5 is not of the same order, at least for larger m = n.

Example 2.11. Consider r = ( p q ) m for p(z) = z, q(z) = z-1 2 with m ≥ 0 an integer. The poles and zeros of r lie in the closed unit disk K = D (or K = [0, 1]) and have a euclidean distance 1 or spherical distance 1/2. However, vec(p m ) 1 = vec(q m ) 1 = 1. In addition we note without proof that

1 (p m , q m ) = 1 (p, q) m = 3 -m and ρ K ( p q ) = ρ {1/3} ( p q ) = 9 4 .
Since ρ K (r) ≤ 2m • ρ K ( p q ) by our previous remark, we can then compare the spherical derivative by

ρ K (r) ≤ 2m ρ K ( p q ) = 2m ρ {1/3} ( p q ) = 9m 2 .
The previous example also shows that the product between spherical derivative and the numerical measure of coprimeness can be very close to 0. Our next statement shows that this product is nicely bounded above. Theorem 2.12.

Let K ⊆ C and r = p q ∈ C m,n [z]. If K ⊆ D or m = n then K,1 (p, q) 2 max ( m vec(p) 1 , n vec(q) 1 ) ≤ 1 ν K (r) ≤ 1 ρ K (r)
.

Theorem 2.12 identifies some particular cases where the bounds of Theorem 2.9(a),(b) are sharper than the bound (2.5) of Theorem 2.5. In addition, in the case of a simple pole z of r = p/q, Theorem 2.12 combined with Corollary 2.8 also implies that p and q numerically relatively prime implies no small residual at z.

Distance of rational functions

Numerical analysis in C m,n [z] requires one to measure distances between rational functions r = p q ∈ C m,n [z] and r = p q ∈ C m,n [z]. As mentioned in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF] several choices are possible. If one is interested in values, then the choice χ K (r, r) := sup z∈K χ(r(z), r(z)) could be the most suitable since the chordal metric measures the euclidean distance of points in the Riemann sphere. On the other hand, if one prefers to define a distance in terms of the coefficients of numerators and denominators, then one should take care of the fact that coefficient vectors are only unique up to a scaling with a complex factor, that is, the norm and the phase. In [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF] the authors made the choice

d(r, r) = min 1 (p, q) 2 vec(p) vec(q) - a ( p, q) 2 vec( p) vec( q) : a ∈ C, |a| = 1 ,
and hence

d(r, r) = min 1 (p, q) 2 vec(p) vec(q) ± 1 ( p, q) 2 vec( p) vec( q)
in the case r, r ∈ R m,n (z) of real coefficients. In [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Theorem 4.1] it is shown that when r is nondegenerate we have

(m + n + 1) -3/2 √
2 cond(S (1) (p, q)) ≤ χ D (r, r) d(r, r) ≤ 2(m + n + 1) cond(S (1) (p, q)).

(2.11) Thus, roughly speaking, the two distances are comparable provided that cond(S (1) (p, q)) is modest, or, in other words, the rational function r is well conditioned.

In view of the preceding statements, it is then natural to wonder whether similar inequalities are kept if we replace cond(S (1) (p, q)) in (2.11) by 1 K,1 (p,q) . The following theorem shows that this is possible for the right-hand side of (2.11). Theorem 2.13. If K ⊆ D or m = n then for any r = p q , r = p q ∈ C m,n [z] we have

K,1 (p, q) • χ K (r, r) ≤ √ 2 (p -p, q -q) 1 .
On the other hand the corresponding statement does not hold for the left-hand side of (2.11) as long as we want to keep constants which are only polynomially growing in m + n. That is, there is not a quantity C(m + n) of order a small power of m + n for which χ D (r, r) ≥ C(m + n) (p -p, q -q) 1 D,1 (p, q), (2.12) Indeed, in the following example we let m = n and show that for each m there are polynomials p m , q m , p m , q m ∈ R m [z], with corresponding rational functions r m and r m , such that, up to some constants,

χ D (r m , r m ) (p m -p m , q m -q m ) 1 ∼ √ m( 3 8 ) m D,1 (p m , q m ) ∼ 9 m cond(S (1) (p, q)) ∼ √ m 27 8 m D,1 (p m , q m ), ( 2
(z) = z m , q m (z) = z-1
2 m be the polynomials from Example 2.11, and suppose these are perturbed as

( p m , q m ) = (p m -ηu m , q m -ηv m )
where

u m , v m ∈ R m-1 [z] are such that q m (z)u m (z) -p m (z)v m (z) = 1
, and η is a small parameter which we will fix later. Then

S (1) (p m , q m ) 1 = vec(p m ) 1 = vec(q m ) 1 = (p m , q m ) 1 = 1. (2.15)
In [3, Examples 5.2 and 5.3] explicit formulas for u m , v m and S (0) (p m , q m ) -1 were derived, allowing the authors to deduce that

vec(v m ) 1 ≤ vec(u m ) 1 ∼ 2 3m-1 √ πm and S (0) (p m , q m ) -1 1 vec(u m ) 1 ∈ [1, 2].
As η is still not fixed, we can now specify that

2 (p -p, q -q) 1 = 2 |η| vec(u m ) 1 = D,1 (p m , q m ), (2.16)
from which, using Example 2.11 and Theorem 2.6(b), we then have for all z ∈ D,

1 (p m , q m ) = D,1 (p m , q m ) = 3 -m with {z},1 ( p m , q m ) {z},1 (p m , q m ) ∈ 1 2 , 3 2 .
Since (p m qm -pm q m )(z) = η, we have

χ D (r m , r m ) = sup z∈D |η| (|p m (z)| 2 + |q m (z)| 2 )(| p m (z)| 2 + | q m (z)| 2 ) ≥ sup z∈D |η| 2 {z},1 ( p m , q m ) {z},1 (p m , q m ) = |η| 3 D,1 (p m , q m ) 2 . Similarly χ D (r m , r m ) ≤ sup z∈D |η| {z},1 ( p m , q m ) {z},1 (p m , q m ) ≤ 2|η| D,1 (p m , q m ) 2 .
Thus, up to some constants,

χ D (r m , r m ) (p m -p m , q m -q m ) 1 ∼ 1 D,1 (p m , q m ) 2 vec(u m ) 1
, which shows our perturbation satisfies equation (2.13).

Finally, notice that (2.16) and (2.15) also imply that

( p, q) 1 ∈ [1 -|η| vec(u m ) 1 , 1 + |η| vec(u m ) 1 ] ⊆ [ 1 2 , 3 2 ],
and thus

d(r m , r m ) ≤ 2|η| vec(u m ) 1 1 -|η| vec(u m ) 1 ≤ 4 |η| vec(u m ) 1 .
Combining, we find that

χ D (r m , r m ) d(r m , r m ) ≥ 1 12 vec(u m ) 1 D,1 (p m , q m ) 2 ,
from which the remaining growth statement (2.14) follows.

Proofs of Theorems

Proof of Theorem 2.3. Denote S = S (0) (p, q) and S * , the matrix obtained by a permutation P of columns of S ( ) (p, q) such that

S * = S ( ) (p, q) • P = S E 0 V , where E ∈ C (n+m)× and V ∈ C (2 )× .
From non degeneracy we know that S is regular and that S ( ) (p, q) is of maximum rank n + m + . Its kernel is therefore of dimension . Let Y ∈ C (n+m+2 )× be a matrix whose columns generate Ker(S * ). We know that the orthogonal projector onto Ker(S * ) is Y (Y * Y ) -1 Y * = Y Y † and so, as S † * S * is the projector onto Ker(S * ) ⊥ we get

S † * S * = I -Y Y † . Let B 2 ∈ C (n+m+2 )× be a matrix such that S ( ) (p, q) • B 2 = 0 I , I identity of order .
We can construct the columns of B 2 in the following way. Let w be the last column of S -1 , that is, Sw = e n+m , where e n+m is the (n + m)th unit canonical vector. In polynomial language, w contains the coefficients of two polynomials u and v of degree n-1 and m-1, respectively, and satisfying the Bezout equation

p(z)u(z) + q(z)v(z) = z m+n-1 .
The columns of B 2 contain the coefficients of z i • u(z) and z i • v(z), i = 1, • • • , . More precisely, the columns have the form

(0, • • • , 0 i-1 , vec(u) T , 0, • • • , 0 -1 , vec(v) T , 0, • • • , 0 -i ) i = 1, • • • .
We set B 1 = P -1 • B 2 and so

B = S -1 0 B 1 (3.1)
is a right inverse of S * since S * B = I. Then

S † * = I -Y Y † B and S † * 2 ≤ B 2 . (3.2)
Let us now bound B 2 . From (3.1) and noting the fact that B 1 2 is bounded by its Frobenius norm of B 1 and remembering that the columns of B 1 are constructed from Sw = e n+m we get

B 2 ≤ S -1 2 + B 1 2 ≤ S -1 2 + √ w 2 ≤ (1 + √ ) S -1 2 .
Thus the second inequality of (2.3) follows.

Let us now prove the first inequality of (2.3). Here we make use of the formula

S E 0 V † = S -1 -S -1 EV † 0 V † .
for the pseudo-inverse of S * . To see this consider the QR decomposition of S and

V S = Q S R S , V = Q V R V 0 = Q V R V ,
where Q S , Q V are unitary matrices and R S and R V invertible since both S and V have maximum rank. Then S * can be written as the product

S * = Q S 0 0 Q V R S E 0 R V , E = Q -1 S E
where the first matrix in the product has linearly independent columns and the second one has linearly independent rows. We then have the pseudo-inverse of S * as

S † * = R -1 S -R -1 S ER -1 V 0 R -1 V Q * S 0 0 Q V = S -1 -S -1 EV † 0 V † .
This trivially gives S -1 2 ≤ S † * 2 and the result follows. Proof of Theorem 2.5. Assume first that p/q is nondegenerate. If we let p, q denote the polynomials with reversed coefficients of p and q, respectively, then we get K,s (p, q) = 1/K,s (p, q). Without loss of generality we may thus suppose min(|z p | , |z q |) ≤ 1. We can write

|p(z q )| = |p(z q ) -p(z p )| ≤ m k=1 |p k | z k q -z k p = |z q -z p | m k=1 |p k | k-1 j=0 z j q z k-j-1 p . (3.3)
Let us suppose |z q | ≥ |z p | and so |z p | ≤ 1. Then by twice applying the Cauchy-Schwarz inequality we obtain

|p(z q )| ≤ |z q -z p | vec(p) 2 m • m-1 j=0 |z q | 2j 1/2 vec(p) 1 m • max(1, |z q | m-1 )
.

The definition of the chordal metric implies |z p -z q | ≤ √ 2 χ(z p , z q ) 1 + |z q | 2 , and using

(1 + |z q | 2 ) m-1 j=0 |z q | 2j = m-1 j=0 |z q | 2j + m j=1 |z q | 2j ≤ √ 2 m j=0 |z q | 2j
we obtain

|p(z q )| ≤ 2 χ(z p , z q ) × vec(p) 2 m • m j=0 |z q | 2j 1/2 vec(p) 1 m • max(1, |z q | m ) .
Thus by the definition of K,s (p, q) K,s (p, q) ≤ 2 χ(z p , z q ) m vec(p) s .

Similarly, if |z p | ≥ |z q | then we get K,s (p, q) ≤ 2 χ(z p , z q ) n vec(q) s and the result then follows.

Finally we remark that the result follows trivially if p/q is degenerate.

Proof of Theorem 2.6. For part (a) we use a Neumann series argument similar to the proof of [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]Lemma 5.1]. Define

E = S (1) (p, q) † S (1) (p, q) -S (1) ( p, q) .
Then by assumption and by comparison with the Frobenius norm • F we have

E 2 ≤ S (1) (p -p, q -q) F S (1) (p, q) † 2 ≤ √ m + n + 1 (p -p, q -q) 2 S (1) (p, q) † 2 ≤ 1/3.
Since also by assumption S (1) (p, q) has full row rank, we find that S (1) ( p, q) = S (1) (p, q)(I -E)

and thus (I -E) -1 S (1) (p, q) † is a right inverse of S (1) ( p, q). This implies that cond(S (1) 1) (p, q)) ≤ 2 cond(S (1) (p, q)).

( p, q)) ≤ S (1) ( p, q) 2 (I -E) -1 S (1) (p, q) † 2 ≤ 1 + E 2 1 -E 2 cond(S (
The other inequality in part (a) follows by symmetry. For part (b), we first notice that, for any z ∈ K,

| {z},s (p, q) -{z},s ( p, q)| = p(z) (1, z, ..., z m ) s , q(z) (1, z, ..., z n ) s s - p(z) (1, z, ..., z m ) s , q(z) (1, z, ..., z n ) s s ≤ p(z) -p(z) (1, z, ..., z m ) s , q(z) -q(z) (1, z, ..., z n ) s s ≤ ( vec(p -p) s , vec(q -q) s ) s = (p -p, q -q) s ≤ K,s (p, q) 2 ≤ {z},s (p, q)
where in the second last inequality we used our hypothesis. Thus 1 2 {z},s (p, q) ≤ {z},s ( p, q) ≤ 3 2 {z},s (p, q) for all z ∈ K. The claim follows by taking the infimum for z ∈ K.

Proof of Theorem 2.9(a),(b). We start by observing that

lim x→z χ(r(z), r(x)) |z -x| = ρ(r)(z), lim x→z χ(r(z), r(x)) χ(z, x) = ν(r)(z),
and thus the suprema in (2.9), (2.10) are bigger than or equal to ρ K (r) and ν K (r), respectively. Following [START_REF] Schiff | Normal Families[END_REF], the spherical metric σ(z 1 , z 2 ) is given by the length of the shortest path on the Riemann sphere joining z 1 and z 2 , and thus

σ(z 1 , z 2 ) = min Γ Γ |dz| 1 + |z| 2
where Γ is any differentiable curve joining z 1 to z 2 , with the minimum being obtained for Γ being the preimage of the shorter of the two circular arcs of radius 1/2 on the Riemann sphere joining z 1 and z 2 . By elementary trigonometry, we can link the chordal metric χ to the spherical metric via

χ(z 1 , z 2 ) ≤ σ(z 1 , z 2 ) ≤ π 2 χ(z 1 , z 2 ). (3.4)
Thus

χ(r(z 1 ), r(z 2 )) ≤ min Γ r(Γ) |dw| 1 + |w| 2 = min Γ Γ |r (z)| 1 + |r(z)| 2 (1 + |z| 2 ) |dz| 1 + |z| 2 ≤ ρ K (r) |z 1 -z 2 | ν K (r) σ(z 1 , z 2 ) .
Combined with (3.4), we obtain the remaining inequalities for establishing (2.9), (2.10).

For the last statements in both (a) and (b) we estimate the distance between pole and zero, with 1 = χ(0, ∞) = χ(r(z p ), r(z q )) ≤ ρ K (r) |z p -z q |, and similarly for ν K (r). Here we have used our assumption that z p , z q ∈ K.

Proof of Theorem 2.12. The inequality ρ K (r) ≤ ν K (r) is an immediate consequence of the definition. Let ν K (r) = ν(r)( z) for some z ∈ K. We consider first the case K ⊆ D and thus | z| ≤ 1. Then ν K (r) ≤ 2ρ(r)( z), where Thus when | z| ≤ 1, the assertion of Theorem 2.12 follows.

If | z| > 1 and thus m = n by hypothesis, we consider as before the reversed polynomials

p(z) = z m p(1/z) ∈ C m [z], q(z) = z n q(1/z) ∈ C n [z],
and observe that r(z) := p(z)/q(z) = r(1/z). Thus ν K (r) = ν 1/K (r) = ν(r)(1/ z) ≤ 2ρ(r)(1/ z), and the assertion follows making use of K,1 (p, q) = 1/K,1 (p, q).

Proof of Theorem 2.13. Using the Cauchy-Schwarz inequality we get (1, z, • • • , z m+n )S (1) (p, q) 1 ≥ 1 |z n | K,1 (p, q) z n+n = K,1 (p, q) |z n | .

|p(z)| 2 + |q(z)| 2 • χ(r(z), r(z)) = |((p -p) q -(q -q) p)(z)| | p(z)| 2 + | q(z) 2 | ≤ |(p -p)(z)| 2 + |(q -q)(z)| 2 ≤ p -p 2 1 max(1, |z| m ) 2 + q -q 2 1 max(1, |z| n ) 2 ≤ √ 2 (p -p, q -q) 1 × 1 if |z| ≤ 1 |z| n if
Then χ(r(z), r(z)) ≤ √ 2 (p -p, q -q) 1 |z| n K,1 (p, q) |z| n = √ 2 (p -p, q -q) 1 K,1 (p, q)

and the result follows.

Conclusions and topics for future research

In this paper we have considered the problem of working with rational functions in a numeric environment, with the particular goal of monitoring the absence of Froissart doublets. Three different parameters were studied, including the euclidean condition number of underlying Sylvester-type matrices depending on some integer , a parameter for determing coprimeness of two numerical polynomials and bounds on the spherical derivative. Each case these parameters sharpen those found in [START_REF] Beckermann | Algebraic properties of robust Padé approximants[END_REF]. Future plans include using our three parameters as penalties for computing rational approximants with the goal of removing such undesirable features. In addition, there are a number of open questions that fall out of our work. All our results are given using a monomial basis. It is of interest to see what can be said in the case of other polynomial bases, for example those based on orthogonal polynomials. It is also natural to look for an interpretation of the quantity ν(r) in terms of residues.

Figure 1 :

 1 Figure 1: Link between the three indicators (Sylvester condition number, coprimeness, bounded spherical derivative) and Froissart doublets of rational functions, together with the numbers of the corresponding statements.

Definition 2 . 1 .

 21 Given an integer ≥ 0 and polynomials p

ρ 2 ≤

 2 (r)( z) = |p ( z)q( z) -p( z)q ( z)| |p( z)| 2 + |q( z)| (|p ( z)| + |q ( z)|) max(|p( z)| , |q( z)|) max(|p( z)| 2 , |q( z)| 2 ) ≤ m p 1 max(1, | z| m ) + n q 1 max(1, | z| n ) max(|p( z)| , |q( z)|) ≤ m p 1 max(|p( z)| , |q( z)|) + n q 1 max(|p( z)| , |q( z)|).

1 ≥

 1 |z| ≥ 1 and m = n .We also have for |z| ≤ 1 that|p(z)| 2 + |q(z)| 2 ≥ max(|p(z)| , |q(z)|) = (1, z, • • • , z m+n )S (1) (p, q) K,1 (p, q) (1, z, • • • , z m+n ) 1 = K,1 (p, q).Using the definition of χ K (r, r) the result follows.For |z| ≥ 1 and m = n we get|p(z)| 2 + |q(z)| 2 ≥ 1 |z n | max(|z n p(z)| , |z n q(z)| = 1 |z n |

(p, q) := inf (p -p, q -q) 1 : ( p, q) ∈ C m [z] × C n [z] have a common root = inf S ( ) (p, q) -S 1 : S a Sylvester-type matrix of not full row rank .

This means that, with z 1 , z

∈ K, also the preimage of the shortest path from z 1 to z 2 on the Riemann sphere belongs to K. Notice that disks and half-planes are spherically convex.


Samuel Bernard

Fabien Crauste

Olivier Gandrillon

Carole Knibbe

David Parsons

Simuscale: A Modular Framework for Multiscale Single-Cell Modelling

Keywords: modelling, multiscale, individual-based model, population dynamics, cell, cell-cell interaction Software architecture 4 3 Plugins 4 4 Numerical

Simuscale is a multiscale, individual-based modelling platform for performing numerical simulations of heterogeneous populations of individual cells evolving in time and interacting physically and biochemically. Models are described at two levels: cellular level and population level. The cellular level describes the dynamics of single cells, as defined by the modeller. Cells have an internal state that includes default properties such as cell size and position, and may also include any other cell-specific state, such as gene or protein expression. The population level describes the mechanical constraints and biochemical interactions between cells. Cells evolve in bounded 3D domain, and can divide or die. Simuscale implements the physical simulator that manages the simulations at the population level. It delegates the details of cellular dynamics to each cell. This makes Simuscale modular, as it can accommodate any number of cell models with the same simulation, including models with different modelling formalisms. Biochemical interactions occur between cells that are in contact with each other, through intercellular signals. Intercellular signals can be known to all or to a subset of the cells only. Simuscale expects an input file describing the initial cell population and numerical options, it runs a simulation over a specified time interval, updating the cell population at given time steps, and generates an output file containing the state of each cell at each time step, and the tree of cell divisions and deaths.

Simuscale: un outil pour la simulation multi-échelles de populations cellulaires

Résumé :

Simuscale est une plateforme de modélisation multi-échelles et individu-centrée pour réaliser des simulations numériques de populations hétérogènes de cellules individuelles évoluant en temps et interagissant physiquement et biochimiquement. Les modèles sont décrits à deux niveaux : cellulaire et populationnel. Le niveau cellulaire décrit la dynamique de la cellule individuelle. Les cellules possèdent un état interne qui comprend notamment des propriétés par défaut comme la taille de la cellule et sa position, et peut aussi inclure d'autres états spécifiques à la cellule, comme l'expression de gènes ou de protéines. Le niveau populationnel décrit les contraintes mécaniques et les interactions biochimiques entre les cellules. Les cellules évoluent dans un domaine 3D borné, et peuvent se diviser ou mourir. Simuscale implémente le simulateur physique qui gère les simulations au niveau populationnel. Elle délègue les détails de la dynamique cellulaire à chacune des cellules. Cette délégation des tâches permet de rendre Simuscale modulaire, en ce qu'elle peut prendre en charge différents formalismes de modélisation. Les interactions biochimiques se produisent entre cellules en contact entre elles, au travers de signaux inter-cellulaires. Les signaux inter-cellulaires peuvent être reconnus en tout ou en partie par un sous-groupe de cellules donné. Simuscale prend en entrée un fichier spécifiant les populations cellulaires initiales et des options de simulation numérique. La simulation met à jour les populations cellulaires à pas de temps discrets sur un intervalle de temps fini, et génère en sortie les trajectoires cellulaires de chaque cellule, ainsi que l'arbre des divisions cellulaire et des morts.

Mots-clés : modélisation, multi-échelles, modèle, individu-centré, dynamique des populations, cellule, interaction cellule-cellule

Introduction

Systems biology is an approach that uses statistical and modelling tools to integrate biological knowledge at lower levels, gene, molecule, cells, with the aim to understand the "big picture" at higher levels: cell, tissues, organisms ( [START_REF] Alber | Integrating machine learning and multiscale modeling-perspectives, challenges, and opportunities in the biological, biomedical, and behavioral sciences[END_REF] and Figure 1a). This raises two challenges. The first challenge is the mismatch between the available biological data or knowledge and the process of interest: there are often missing links that prevent the bridging of scales. Multiscale models are computational models that provide an explicit description of the links between the different scales, often with interactions that are bi-directional [START_REF] Bernard | How to Build a Multiscale Model in Biology[END_REF]. Many problems concern the determination of cell fate (differentiation, survival, death, proliferation) from biomolecular data, and the nature of cell decision mechanisms. The second challenge is heterogeneity, when outliers cannot be taken into account by averaged models, but are able to tilt the big picture. The main problem we face in computational implementation of multiscale models is linking or bridging the scales. Assuming that a cell can be described (observed) by a vector of low level markers X in a set Ω, what can be said about the fate of the cell, or of the tissue to which it belongs? Even with a detailed model of the molecular processes going on inside the cell, it is difficult to infer the fate of the cell unambiguously.

In this report we present the simulation platform Simuscale for multiscale and integrative modelling in systems biology. The development of Simuscale1 addresses the two challenges mentioned above by providing a simplified single-cell-based simulation environment (Figure 1b). Simuscale is a versatile multiscale modelling plateform for efficient numerical simulation of large, heterogeneous populations of interacting cells. Cells are represented by oriented visco-elastic balls that evolve under mechanical constraints in a three-dimensional domain W ⊂ R 3 . Each cell possesses a volume, a position and an orientation. Additionally, each cell can be endowed with user-defined intrinsic intracellular dynamics. The intracellular states of each cell can be coupled to other cells in contact via biochemical signals, as if expressed on the surface of the cells. Macroscopic cell fate (division, death, growth, differentiation, or movement) is entirely dictated by the intracellular state of the cell. Physical and biochemical interactions evolve in a 3D environment. Simuscale takes the form of a core numerical simulation engine on top of which biological model plugins can be added by the user. Plugins are meant to be interoperable, so that different plugins (say, one for tumor cells and one for immune cells) can be used together in the same simulation. This feature relies on the fact that all cells share common properties: position, size, orientation, mechanical properties, and biochemical signals. The intracellular dynamics is encapsulated in each cells, so different types of formalism can co-exist: Boolean networks, ordinary or stochastic differential equations, piecewise-deterministic Markov processes. This reduces the need for re-implementing existing models. There are several existing modelling frameworks, such as CompuCell3D [START_REF] Graner | Simulation of biological cell sorting using a twodimensional extended Potts model[END_REF], PhysiCell [START_REF] Ghaffarizadeh | PhysiCell: an open source physics-based cell simulator for 3-D multicellular systems[END_REF], Chaste (Cancer, Heart and Soft Tissue Environment), TiSim/CellSys INRIA/IfADo [START_REF] Hoehme | A cell-based simulation software for multi-cellular systems[END_REF], and Tissue Forge [START_REF] Sego | Tissue forge: Interactive biological and biophysics simulation environment[END_REF]. By contrast, Simuscale is not focused on mechanical interactions. It rather focus on intracellular dynamics and cell-cell interaction, by providing a modular and plug-and-play interface to run simulations with several interacting cell models.

Software architecture

Simuscale numerically simulates the time evolution of a heterogeneous population of interacting cells within a bounded spatial domain. Simulations are performed iteratively in time with discrete time steps. At each time step, intercellular signalling and mechanical interactions are first computed for each cell. The internal state of each cell is then updated, including cell size and position, and cell death and cell division processed. The update is performed synchronously.

Simuscale is structured around a simulation core that performs model initialization, time iteration, and generates the simulation output. Iterations are synchronous: mechanical, biochemical interactions and cell fate (death or division) are first computed, then cells are updated based on computed interactions. Cellular updating is delegated to user-defined plugins that implement the details on how cells should evolve.

Because the focus is not on the mechanical behaviour, the simulation core includes predefined cell movement behaviours. The user can easily implement three types of movement: non-moving cells (immobile), cells that can move due to mechanical constraints, but have no motility of their own (mobile), and cells that can move on their own (motile). Movement rules for motile cells is user-defined.

Biochemical signals are mediated by a signalling interface common to all plugins. The signalling interface consists in a list of signals that can be expressed and detected by each cell. During the core update, each cell scans the subset of signals expressed by neighbouring cells and integrates the total signal strength. There is no limit to how many different signals a cell can express or detect. Each plugin defines the subset of signals they can use, and there is not need for the subsets of different plugins to intersect.

The numerical iterations are performed over a predefined time interval t ∈ [0, T ]. During iterations, a partial set of the state of each cell is recorded in a text file called trajectory.txt, in a tall format, i.e. each row corresponds to a single cell at a given time. Cell state includes: time, a unique cell identifier, spatial position, and radius. Optionally, the orientation and the expression of any subset of signals (expressed or not by the cell) can be recorded. The text file format was chosen so that it can easily be parsed into tables and data frames for post-processing by standard packages in Python or R, such as panda or dyplr. In addition to the trajectory file, a cell lineage tree is recorded in Newick format [START_REF] Felsenstein | PHYLIP (phylogeny inference package), version 3.5 c[END_REF][START_REF]Newick format[END_REF] in the file newick.tree, and in a graph format (TGF) [START_REF]Trivial graph format[END_REF] in the file treelist.tgf. Numerical simulation algorithms are detailed below.

Plugins

Cellular models are implemented in plugins. A plugin is a class derived from the class Cell that implements the internal dynamics, the movement behaviour, intercellular signalling, growth, and cell division and death (Figure 2).

Movement behaviour

Each cell is modelled as a visco-elastic ball with a soft shell and an incompressible core. The soft shell is the intersection of two concentric balls with radii r and s, 0 < r < s. The core is the ball of radius r. The radius r is called the internal radius, and the radius s is call the external radius. The volume of a cell can change with time, but we assumed that the ratio of the internal to the external radius is a user-defined constant parameter R.

Cells have one of the three movement behaviours: immobile, mobile, or motile. Cell with immobile behaviour do not move. This behaviour is suited for cells that play a static role, and do not divide. Cells with mobile behaviour can move when forces are applied on them, but do not generate a force on their own. Cells with motile behaviour can move when forces are applied on them, and can also generate a force on their own, as defined by the model.

To model the mechanical interaction between cells, we use a capped Lennard-Jones type potential in a high-viscosity environment, i.e. cell displacement is proportional to the force applied to it. The force exerted between two cells from the Lennard-Jones potential is slightly repulsive (positive) when cells are close, but not in contact. The force becomes attractive, i.e. takes a negative value, when cells get in contact with each other, and become repulsive again when the distance between the centres of the two cells become less than the sum of the internal radii. Let r 1 and r 2 be the internal radii of the two cells, and δ the distance between the two cells. The force due to physical cell contact is

F = -24 2 σ 12 δ 13 - σ 6 δ 7 ,
with parameters = 0.002, σ = (r 1 + r 2 )/L, and L ∼ 1.12. When the ratio of the internal to the external radii is around 0.9, cells tend to stick together once they become in contact. When the distance between the centres of the two cells is less than the sum of the internal radii, the incompressible cores overlap, and the force quickly becomes repulsive. When δ is small, the Inria force F can become large. To avoid unnecessary numerical instabilities when solving equations of movement, we capped the force F by a constant value F max = 0.5,

F c = cap(F, F max ),
where the capping function defined by

cap(x, M ) = x, |x| ≤ M, sign(x)M, |x| > M.
The mechanical movement equations ẋ = F , x ∈ R 3 , are solved with an explicit updating scheme.

Intercellular signals

Intracellular dynamics

The intracellular dynamics is implemented as a blackbox, and is left to the user to define. The advantage of this approach if that in principle any computable model can be implemented. The drawback is the time cost of implementation.

Cell fate

Cells can divide and die. Rules for death and division are implemented in the plugin. At division, a daughter cell is created from the mother cell. The daughter cell shares the same class as the mother cell. The initial state of the daughter cell depends on the mother cell only. The state of the mother cell can be reset at division to allow the partitioning of biological material, for example. At death, the cell is instantaneously deleted from the simulation.

Cell type

The single-cell approach used in Simuscale means that cells do not have to have pre-defined cell types. However it can be useful to tag cells into a discrete classes, for computational or analytic purposes. The file CellType.values.in contains a list of cell types that can be used. As with the intercellular signal, any number of unique cell types can be added. Cell types must be valid C++ identifiers.

Extracellular environment

Currently, simulations are exclusively cell-based, and the extracellular environment is empty. However, non-cell extracellular particles (such as macro-molecules) can be implemented as a specific plugin.

Numerical simulations

Launching. Plugins are located in the directory src/plugins. As plugins are written directly as source code, a new compilation is necessary after changes in plugins. Simulation parameters are contained in a file named param.in, in which initial population sizes and types can be set. Simulation is launched with the executable build/bin/simuscale. Numerical simulations are performed iteratively with a discrete time step parameter DT hours, according to the flowchart shown in Figure 3.

Start. At the beginning of the simulation, cell populations defined in param.in are created. At cell creation, the constructor defined in the plugin is called. By default, spatial initial positions are chosen randomly. User-defined initial positions can be set in two ways: 1) by overriding the initial position in the constructor, and 2) by providing the positions in a text file, detailed below.

Initial cell populations are specified by the keyword ADD_POPULATION, followed by six mandatory parameters: number of cells, cell type, plugin, movement behaviour, doubling time, minimal volume. Optionally, the name of a file containing the initial positions and radii of each initial cell can optionally be provided. The initial position file is a space separated file with six columns:

L(int) X(float) Y(float) Z(float) R(float) C(int)
Only columns 2 to 5 (initial position (x, y, z) and the initial external radius c) are used at the moment. Cells are initialized in row order. The total number of rows must be at least equal to the number of cells in the population being added.

Iteration. Iterations are synchronous: mechanical forces and intercellular signal from neighbouring cells are are first computed. Cells are then updated independently: movement, orientation, internal status including the decision for cell death and cell division, followed by cell death and cell division. Cell death has priority on cell division. If a cell is marked for death and division, division will not occur.

End. Simulation ends when t ≥ T . The complete cell lineage tree is written in two formats, newick (newick.tree) and TGF (treelist.tgf). The trajectory file can be visualised with the companion tool build/bin/view.

Parameter file

The parameter file param.in is necessary to run the simulations. The general syntax is Inria KEYWORD VALUE1 VALUE2 ... Lines starting with # and empty lines are ignored. Keywords are listed in Table 1. Any number of lines with the keywords ADD_POPULATION and SIGNAL can be present in the parameter file.

Output files

The main output file is the text file trajectory.txt. It contains the trajectory of each indiviudal cell over the interval t ∈ [0, T ]. Lineage trees are stored in newick.tree and treelist.tgf. The file normalization.txt is a space-separated value file that contains the extremal values of each signal present in the trajectory file the entire time interval: each line has the form signal max min. The normalization file is intended to be used with the companion script view, which can display a 4D rendering (3 spatial dimensions plus time) of the trajectory file. Each cell is identified by a unique cell ID, starting at ID=1, that remains the same across its life span. At division, a new cell is created with a unique ID.

The trajectory file is text file beginning with a header, and followed by a space-separated value data section. The header is formed by the lines starting with one of the characters !, # or $. The data section starts with the first line starting with a floating-point number. The line starting with a ! specifies the world size in the x, y, z dimensions. The lines starting with a # or $ describe each column in the data section. Column description has the form column number : description. Lines starting with the character # describe columns that have a fixed meanings: time, number of cells in the population, cell ID, x position, y position, z position, x orientation, y orientation, z orientation, external radius. Orientation are written to the trajectory file only if WRITEORIENTATION is set to 1. Lines starting with $ describe signal column. Only signal specified in the parameter files are included in the trajectory file, even though any signal can be used in the simulation. The data section has a fixed number of columns as specified in the header. Each row records the state of a single cell at a given time (it is a tall table format with respect to time).

The newick file is a text-based file encoding a binary tree that records the whole cell lineage. Unlike many graph formats, the newick format provides a way to specify the length of each edge with the syntax node:length. Inner nodes (non-leaf, non-root nodes) represent cell division, with one branch for the mother cell and one branch for the daughter cell. Leaf nodes represent the end of existence of a cell (either death or end of simulation). A cell division is recorded as (descendants)cellID>motherID:ageOfMother where descendant is the tree with a root at cellID. Cell deaths are recorced as tree leaves cellID:ageAtDeath. Cells that are alive at the end of the simulation are recorded the same way, with cellID:ageAtT. A tree is a comma separated list of trees enclosed by parentheses and followed by cellID>motherID:ageOfMother. Because cells initially present at t = 0 do not have ancestors, a "root" cell with ID 0 is added to the tree to make it rooted. Therefore the newick tree always ends with ...>root:0.00)2>root:0.00)1>root:0.00;

The annotation is made compatible with the phylogenetic plot method from the ape package in R [START_REF]Analyses of phylogenetics and evolution[END_REF].

The trivial graph format is a text-based adjacency list format, i.e, the graph is represented as a list of edges. The file is divided in two parts: node definition and edge definition. Node definitions and edge definitions are separated by a line containing a single #. In the first part, each line defines a node: node(uint32) label(string without space). The node is the cell ID, and the label is <cell_type>-ID. In the second part, each line defines a directed edge as a pair of nodes. The TFG format does not provide a native way to specify edge lengths, but 5 Test case: circadian clock-regulated tumour growth

As a test case, we re-implemented a model for the control of the cell cycle by the circadian clock [START_REF] El Cheikh | A multiscale modelling approach for the regulation of the cell cycle by the circadian clock[END_REF]. In mammals and other organisms, cell division is known to occur at preferential times of the day. The circadian clock is an autonomous biological rhythm that has been characterized across all life domains. The circadian clock is generated and maintained by gene regulatory loops that modulate the activity of transcription factors and their downstream targets. The cell cycle genes p21 (controlling cell growth) and Wee1 (controlling mitosis) are under direct control of the circadian clock [START_REF] Farshadi | Molecular links between the circadian clock and the cell cycle[END_REF]. The circadian clock and the cell cycle models, implemented as coupled systems of ODEs, were adapted from previously published models [START_REF] El Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF][START_REF] El Cheikh | A multiscale modelling approach for the regulation of the cell cycle by the circadian clock[END_REF]. The circadian clock module is described by seven variables (Figure 4, circadian clock module), and the cell cycle by two variables (Figure 4, cell cycle module). The circadian clock transcription factor BMAL1 acts directly on the cell cycle variable X, which has a role similar to Wee1. BMAL1 activates X, which in turn blocks the transition to mitosis and cell division. As the activity of BMAL1 decrease during the day, X becomes less inhibited, and cell division can proceed. As an additional mechanism, the circadian clock synchronizes neighbouring cells through activation of the Per2 and Cry circadian clock genes via a Clock signal. As a result, cells tend to divide at preferential times. Cell death occurs when local cell density rise, limiting the total population (Figure 4, cell The initial configuration consisted in 17 cancer cells located near the center of the computational domain (the World). An 18th cell, a killer cell, was added. Killer cells move randomly. Upon contact with cancer cells, the killer cell reduces its motility, increase the probability of clock cells to die, and starts proliferating. The numerical simulation ran for t in [0, 480] hours (Table 2). At t = 480h, the total cell number was 23,706 (cancer cells: 21,202; killer cells: 2504). The Clock and Cycle signals maintained a moderate amount of coherence during the first six days (144 hours), but Clock synchrony (and consequently Cycle synchrony) disappeared once the killer cell made contact with cancer cells and started proliferating (Figure 5a,b and Figure 6a,b,c). The Death signal was restricted to the zone in contact with the killer cells (Figure 6d).

Conclusion

Preliminary tests show that Simuscale is viable for modelling cell population in realistic settings. While the current release is functional, there is a number of planned features to be added in future releases. First, there is a need for a better user-interface, with a clearer separation between the core simulator and the plugins. Moreover, an integrated initial configuration tool could be helpful. Currently, modifying the initial state of the population requires changes in several files. Second, although there are many advantages to using a single-cell modelling approach, the lack of an extracellular environment can limit certain applications. For example, long-range intercellular interaction is not possible, because only cells in contact with each other can communicate. Adding a Gaussian field to cells could mimick the diffusion of molecular signals at longer ranges, at a light computational cost, thank to fast Gauss transform algorithms [START_REF] Greengard | The fast Gauss transform[END_REF]. Numerical optimisation could take several form. Mixed-precision numerical schemes exploit reduced precision arithmetics to speed up computation without reduction in numerical accuracy2 . This could be particularly relevant in the context of parallelisation on GPUs, where large amount of data transfer can severely limit speed [START_REF] Gregg | Where is the data? Why you cannot debate CPU vs. GPU performance without the answer[END_REF]. Working with low precision arithmetics, such as half precision, could have the advantage of benefiting from the hardware acceleration on modern GPUs developed for machine learning. Finally, developing concurrent or post-mortem (once the simulation ended) data analysis tools will be useful for model or parameter inference and optimisation. Tools for macroscopic/continuum approximation could also be valuable to gain analytical insight, or to reduce the numerical complexity of the models [START_REF] Matthew J Simpson | Reliable and efficient parameter estimation using approximate continuum limit descriptions of stochastic models[END_REF].

Figure 1 :Figure 2 :

 12 Figure 1: Scales (a) and modelling approaches (b) in systems biology. Simuscale is centrered around the cell (grey area in a) and upper-left corner in b)).

Figure 3 :

 3 Figure 3: Flowchart of the numerical simulation in Simuscale .

Figure 4 :

 4 Figure 4: Circadian clock-driven cell cycle model. Adapted from reference [11].

Figure 5 :

 5 Figure 5: Tumour growth under a circadian clock control. a) Clock signal. The Clock signal synchronizes cells together, and drive the cell cycle. Clock synchronization decreases with tumour growth. b) Cycle signal. High levels indicates that the cells are more likely to divide. c) Death signal. Cell death occurs when the Death signal is above 1.0. States at time points 0, 72, 144, 216, 288, 360, 432 and 480 h. Cell types (Cancer and Killer) are best distinguished in panel c), where cancer killer cells are in dark blue and cancer cells have colours corresponding to positive Death expression, from light blue to red.

Figure 6 :

 6 Figure 6: Total cell number a), mean signal expression levels b, Single-cell Clock signal expression levels for the initial 400 cells, and the spatial distribution of the cell cycling status d. In (d, dying cancer cell (red) co-localize with the killer cells (blue).

  Cells interact by expressing intercellular signals to neighbouring cells. The intercellular signals expressed by a cell is a subset of a common set of intercellular signals. Neighbouring cells can pick up these intercellular signals. There is no limit to the number of intercellular signals. The full list of intercellular signals is maintained in a master file InterCellSignal.values.in. Signals can be defined in the list as needed. A signal name must be a valid C++ identifier.

Table 1 :

 1 Description of keywords in the parameter file param.in.

	RT n°0520

Table 2 :

 2 Simulation files for the test case. density module). For the purpose of the test case, the cells endowed with this model are called cancer cells. The file param.in is

	File		Description
	currentdir/param.in	Parameter configuration file
	src/plugins/Cell_SyncClock.cpp clock cell implementation
	src/plugins/Cell_SyncClock.h	header file for clock cells
	src/plugins/Cell_Killer.cpp	killer cell implementation
	src/plugins/Cell_Killer.h	header file for killer cells
	PRNG_SEED	1421746927
	MAXTIME	480.0
	DT	0.1
	ADD_POPULATION	17 CANCER SYNC_CLOCK MOTILE 40 1.0 init_clock_4.txt
	ADD_POPULATION	1 KILLER KILLER MOTILE 8 1.0 init_killer.txt
	USECONTACTAREA	1
	WRITEORIENTATION	1
	R_RATIO	0.90
	SIGNAL CLOCK
	SIGNAL CYCLE
	SIGNAL DEATH
	SIGNAL KILL

Simuscale is hosted at https://gitlab.inria.fr/bernard1/simuscale.RT n°0520

RT n°0520

Inria exploratory project ExOdeInria

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-0803