Nidhi Kaihnsa

Yue Ren

Mohab Safey

El Din

AND Johannes W R Martini

Cooperativity, absolute interaction, and algebraic optimization

published or not. The documents may come

Introduction

Interaction between components is a fundamental feature of biological systems. While a simple system of independent subunits is completely defined by its subunits, a complex system with interactions is more than the sum of its parts. A classical example of a small biological system with non-trivial interaction is hemoglobin with its four binding sites for oxygen [BHK04; [START_REF] Barcroft | The combinations of haemoglobin with oxygen and with carbon monoxide. II[END_REF][START_REF] Hill | The combinations of haemoglobin with oxygen and with carbon monoxide. I[END_REF]. The ligand oxygen binds to the four binding sites of the (target) molecule hemoglobin and the interaction can be seen on the overall (isotherm) binding curve, which relates the average ligand saturation inside the system to the ligand concentration outside the system at constant temperature [START_REF] Hill | Cooperativity theory in biochemistry: steady-state and equilibrium systems[END_REF]. While the binding curve of an independent system is linear, the binding curve of hemoglobin has a sigmoidal shape [START_REF] Barcroft | The combinations of haemoglobin with oxygen and with carbon monoxide. II[END_REF][START_REF] Hill | The combinations of haemoglobin with oxygen and with carbon monoxide. I[END_REF], see Figure 1. The sigmoidal shape implies that the extremal states of full and zero saturation are more stable than the intermediate states of partial saturation. This phenomenon is commonly referred to as cooperativity, named after the intuition that bound ligands affect, either positively or negatively, the chances of new ligands binding to the still open sites (see [START_REF] Stefan | Cooperative binding[END_REF] for a review and interpretations). Cooperativity is ubiquitous in nature. It is an essential trait for transport molecules, and its importance in the formation of multi-protein complexes [START_REF] Roy | Cooperative binding mitigates the high-dose hook effect[END_REF], in general signal transduction processes [Sal+16; LFSR09; MSS16] and in the regulation of noise [GMD09; MMD13] has been discussed.

When investigating cooperativity, one interesting problem is to connect macroscopic behavior to microscopic behavior. Here, macroscopic refers to the measured average saturation at varying ligand concentration, and microscopic refers to which individual sites are saturated and how the different sites interact. An example is the aforementioned interpretation of occupied sites increasing the probability of ligands binding to unoccupied sites [START_REF] Stefan | Cooperative binding[END_REF][START_REF] Onufriev | Decomposing complex cooperative ligand binding into simple components: connections between microscopic and macroscopic models[END_REF]. Ideally, such connections between microscopic and macroscopic behavior provide insight on the potential structure and the inner workings of the molecule.

However, the overall binding curve alone does not allow for much inference on the individual site binding curves, similar to how a large sum allows for little inference on its summands. Given a fixed overall binding curve, there generally is no unique molecule, defined by binding and interaction energies, with the prescribed macroscopic binding behavior, but infinitely many. Thus many works feature additional constraints which reduce the ambiguity, such as considering only small system with one or two binding sites [START_REF] Hunter | What is cooperativity?[END_REF] or symmetric systems with identical binding sites [START_REF] Abeliovich | On Hill coefficients and subunit interaction energies[END_REF][START_REF] Rong | Homogeneous molecular systems are positively cooperative, but charged molecular systems are negatively cooperative[END_REF]. These constraints have the added benefit of simplifying the mathematical problem, allowing for easier inference on interaction energy terms, conditional probabilities, correlations or other quantities. Ideally, constraints should be motivated by knowledge on the chemical structure of the molecule.

A biologically and mathematically interesting question is for instance whether there exist a molecule with independent sites for a fixed overall binding curve. Assuming evolutionary pressure on the transport properties of hemoglobin, it should not be possible to model the system's overall binding curve through an independent system, since otherwise the tetrameric structure of hemoglobin would hold no advantage over four monomers.

The decoupled sites representation [OCU01; MU13] is a central theorem in this context, which states that every overall binding curve can be generated by a unique hypothetical molecule with possibly non-real binding energies and independent binding sites. The binding energies are determined by the roots of the binding polynomial of the overall binding curve. In particular, the realness of the roots show whether interaction is necessary to achieve the overall binding behavior and thus provide a qualitative definition of cooperativity.

In case interaction is required to generate the overall binding curve, a natural question is how to quantify the degree of cooperativity. This is classically done using the maximal slope of the Hill plot, which is the plot of binding curve on a logarithmic x-axis (see Figure 6). A new approach [START_REF] Martini | A measure to quantify the degree of cooperativity in overall titration curves[END_REF] is to look at the minimal interaction required to generate the considered overall binding curve. The minimal interaction quantifies the minimal deviation from an independent system that is necessary to explain the observed overall binding behavior.

The manuscript in hand studies the problem of finding molecules generating a given overall binding curve with minimal interaction. We show that looking for the minimal interaction energy required to generate a given overall titration behavior, leads to an algebraic optimization problem which can be partially tackled using software such as SCIP [START_REF] Gleixner | The SCIP Optimization Suite 6.0[END_REF]. We calculate the minimal molecules of various binding curves in existing datasets [START_REF] Connelly | Analysis of zeros of binding polynomials for tetrameric hemoglobins[END_REF][START_REF] Ikeda-Saito | Thermodynamic properties of oxygen equilibria of dimeric and tetrameric hemoglobins from Scapharca inaequivalvis[END_REF][START_REF] Imai | Analyses of oxygen equilibriums of native and chemically modified human adult hemoglobins on the basis of Adir's stepwise oxygenation theory and the allosteric model of Monod, Wyman, and Changeux[END_REF]. These binding curves describe the oxygen binding of hemoglobins under different chemical treatments and under different temperatures. We compare their minimal interaction with the maximal slope of the Hill plot, and discuss similarities and discrepancies with a view towards the shapes of the binding curves.

Recapitulation of the mathematical framework

In this section, we briefly review the model for ligand binding based on the grand canonical ensemble of statistical mechanics and recall the notion of minimal absolute interaction from [START_REF] Martini | A measure to quantify the degree of cooperativity in overall titration curves[END_REF]. We refer the reader to [Ben13; WG90; Hil85] for detailed expositions on the model. Definition 2.1 A hypothetical molecule W with n sites is a complex point, whose coordinates are indexed by the subsets of [n] := {1, . . . , n} and with w ∅ = 1:

W = (w I ) I⊆[n] ∈ C 2 n .
A (real) molecule is a hypothetical molecule whose coordinates are real and positive.

We refer to w I as binding energies if |I| = 1 and interaction energies if |I| > 1. For the sake of brevity, we abbreviate w {i 1 ,...,ir} to w i 1 ...ir for i 1 < • • • < i r (see Fig. 2).

Note that we refer to w i 1 ...ir as energy even though it is actually representing the exponential of an energy difference [Ben13; WG90; Sch75]. To be precise:

• The binding energy w i encodes the difference between the free energies of the completely unsaturated state and the state with a single ligand bound to site i. • The interaction energy w i 1 i 2 encodes the difference between the free energy when sites i 1 , i 2 are occupied and the sum of the energies when either i 1 or i 2 is occupied. • Analogously, w i 1 ...ir encodes the difference between the free energy when sites i 1 , . . . , i r are occupied and the sum of all energies when a proper subset of i 1 , . . . , i r are occupied.

As such, the only "sensible" values for w i 1 ...ir are positive real numbers. However, mathematically it is easier to work with hypothetical molecules, and they allow for interesting phenomena such as the decoupled sites representation [START_REF] Onufriev | A novel view of pH titration in biomolecules[END_REF][START_REF] Martini | A mathematical view on the decoupled sites representation[END_REF]. In this article, we will mainly focus on real molecules. Definition 2.2 The overall binding curve of a real molecule W = (w I ) I⊆[n] ∈ (R >0 ) 2 n with n sites is defined as The binding polynomial Φ(W) is the univariate polynomial of degree n in its denominator,

Ψ(Λ) = na n Λ n + (n -1)a n-1 Λ n-1 + • • • + a 1 a n Λ n + a n-1 Λ n-1 + • • • + a 1 + a 0 ,
Φ(W) := a n Λ n + • • • + a 1 Λ + a 0 .
The binding polynomial Φ uniquely determines the overall binding curve Ψ.

Example 2.3 Let W = (w I ) I⊆[3] ∈ (R >0 ) 2 3 be a molecule with 3 sites. The binding polynomial Φ(W) = a 3 Λ 3 +a 2 Λ 2 +a 1 Λ+a 0 is a real univariate polynomial of degree 3 whose coefficients are given by

a 0 = w ∅ = 1, a 1 = w 1 + w 2 + w 3 ,
a 2 = w 1 w 2 w 12 + w 1 w 3 w 13 + w 2 w 3 w 23 , a 3 = w 1 w 2 w 3 w 12 w 13 w 23 w 123 .

Note that two different molecules may have the same binding polynomial and thus the same binding curve, e.g.,

w ∅ w 1 w 2 w 3 w 12 w 13 w 23 w 123 W = ( 1 3 1 1 1 1 2 2 ), W = ( 1 2 2 1 1 1 1 3 ).
Hence, the map Φ which maps a molecule with n sites to its binding polynomial is not injective.

Computing hypothetical or real molecules W with a fixed binding polynomial P can be a hard task, though there are well established algorithms and software for both. Let M C (P ) resp. M R (P ) denote the set of hypothetical resp. real molecules W with Φ(W) = P . Then M C (P ) is a algebraic variety, i.e., a set of complex points given by polynomial equations, while M R (P ) is a semialgebraic set, i.e., a set of real points given by polynomial equations and inequalities.

Computing a hypothetical molecule in M C (P ) is generally easier than M R (P ). It can be done purely numerically, e.g., using numerical algebraic geometry and software such as bertini, PHCPack, Hom4PS or juliahomotopycontinuation.

Computing a real molecule in M R (P ) is a significantly more involved task. There are many approaches to it and most require some symbolic computation, such as RAGlib which uses a mix of symbolic and numerical computations.

However, both M C (P ) and M R (P ) are of high dimension and degree [START_REF] Ren | Decoupled molecules with binding polynomials of bidegree (n,2)[END_REF], so that it is unclear what information can be read off a randomly computed molecule. Rather, we need sensible constraints which gives us molecules of interest. Since the aim of this paper is to study cooperativity in terms of the interaction between the binding sites, we introduce the notions of absolute interaction and minimal interaction:

Definition 2.4 The (absolute) interaction of a molecule W = (w I ) ∈ (R >0 ) 2 n is given by W := |I|>1 max w I , w -1 I .
Since w I represents the exponential of a difference of actual energies, which can be positive or negative, max(w I , w -1 I ) represents the exponential of the absolute value of that difference. Hence, the absolute interaction W represents the exponential of the sum of the absolute values of all energy differences. The minimal value W can adopt is 1 which corresponds to a minimal physical interaction energy of 0. Definition 2.5 Given a binding polynomial P of degree n, the minimal interaction is

P := min W W ∈ (R >0 ) 2 n with Φ(W) = P .
We call a molecule W minimal, if W = Φ(W ) .

In other words, the minimal interaction P is the minimal interaction required to generated the binding polynomial P . Thus it quantifies the minimal deviation from an independent system that is required to explain the observed overall binding behavior. The following proposition states that the minimum exists.

Proposition 2.6 ([Mar17, §4]) For any univariate polynomial P of degree n with positive coefficients there exists a molecule W ∈ (R >0 ) 2 n with n sites such that Φ(W) = P and W = P .

Since cooperativity is a property of the binding curve which emerges from the interaction between the sites, the minimal interaction of a binding polynomial is a natural candidate for quantifying cooperativity. In addition, [START_REF] Martini | A measure to quantify the degree of cooperativity in overall titration curves[END_REF]§4] shows that it has several properties that are desirable of a quantifier.

Example 2.7 Consider the following binding polynomials:

P 1 :=4Λ 3 + 3Λ 2 + 2Λ + 1, P 2 :=6Λ 3 + 7Λ 2 + 4Λ + 1 = (2Λ + 1)(3Λ 2 + 2Λ + 1), P 3 :=6Λ 3 + 11Λ 2 + 6Λ + 1 = (Λ + 1)(2Λ + 1)(3Λ + 1).
A brief computation reveals that the minimal absolute interactions are

P 1 = 13.50, P 2 = 3, P 3 = 1.
The computation for P 1 is explicitly shown in Example 3.3. Figure 3 illustrates the minimal molecules for each binding polynomial. Since P 3 factorizes into three real linear factors, its minimal molecule has only trivial interactions.

Computing the minimal interaction

In this section, we consider the problem of computing the minimal absolute interaction for a given binding polynomial P :

minimize W W subject to Φ(W) = P (1)
For the sake of clarity, we will restrict ourselves to the case of molecules with 4 binding sites, though our arguments also apply to the general case. We will discuss what makes Problem (1) challenging, and what can or cannot be done towards solving it from a practical point of view.

3.1. The initial problem. Problem (1) is commonly called the standard form of a continuous optimization problem with W being the objective function and Φ(W) = P the constraints. The set of all molecules W which satisfy the constraints is the feasible set.

Generally speaking, there are many approaches for solving optimization problems depending on the type of the objective function and the constraints, which may also include inequalities. If both are linear, the problem is also known as a linear program. If the objective function is quadratic and the constraints are linear, the problem is also called a quadratic program. Both types of optimization problems are well-studied with plenty of literature and software solutions. If both objective function and constraints are polynomial, the problem is a polynomial optimization problem or an algebraic optimization problem. This is the type of problem that needs to be solved for Problem 1, and solving it generally involves using sums of squares approaches [Las01; Par03] or critical point methods [START_REF] Greuet | Probabilistic algorithm for polynomial optimization over a real algebraic set[END_REF].

For molecules with 4 binding sites specifically, Problem (1) is of the following form for a given binding polynomial P = a 4 Λ 4 + a 3 Λ 3 + a 2 Λ 2 + a 1 Λ: minimize (2) Problem (2) is extremely hard to solve because of two reasons:

(w I ) I⊆[4] |I|>1 max w I , w -1 I subject to a 1 = w 1 + w 2 + w 3 + w 4 , a 2 =
• the objective function contains 11 maxima,

• the constraints contain polynomials of degree up to 15.

In the following, we will briefly discuss the two challenges separately.

3.2.

Resolving the maxima in the objective function. The maxima in the objective function of Problem (2) can be resolved in two ways. One standard way is to reformulate the objective function as a maximum of polynomial functions and linearize the objective function through additional constraints: Abbreviating [4] >1 := {I ⊆ [4] | |I| > 1}, we define for each set of interactions

I ⊆ [4] >1 f I ((w I ) I⊆[4] ) := I∈( [4] >1 ) I∈I w -1 I • I∈( [4] >1 ) I ∈I w I , so that (w I ) I⊆[4] = max I⊆( [4] >1 ) f I ((w I ) I⊆[4] ).
Then Problem (2) can be written as: minimize (

(w I ) I⊆[4] r subject to f I ((w I ) I⊆[4] ) ≤ r for all I ⊆ [4] >1 , a 1 = w 1 + w 2 + w 3 + w 4 , a 2 =
) 3
This resolves all maxima in the objective function at the cost of introducing 2 11 additional constraints, one for each I ⊆ [4] >1 .

Alternatively, one can decompose the space of molecules (R >0 ) 2 4 into 2 11 regions in which either w I ≥ 1 or w I ≤ 1 for all I ⊆ [4] >1 . For each I ⊆ [4] >1 , we define a region

O I := (w I ) I⊆[n] ∈ (R >0 ) 2 4 | for all I ∈ [4] >1 : w I ≥ 1 if I / ∈ I and w I ≤ 1 if I ∈ I so that (R >0 ) 2 4 = I⊆( [4] >1 ) O I and W = f I (W ) for W ∈ O I .
Finding the minimal absolute interaction in O I hence becomes a straight-forward polynomial optimization problem with few constraints.

minimize

W∈O I f I (W) subject to Φ(W) = P (4)
However, this approach requires solving 2 11 polynomial optimization problems, one for each region O I , which is not easily feasible. One can reduce the number of regions that require consideration slightly by exploiting the intrinsic symmetry of the problem.

(x,y)∈R 2 >0 max(x, x -1 ) • max(y, y -1 ) subject to a 1 = x + y, a 2 = x • y.
The space R 2 >0 can be decomposed into 4 regions on which the objective function is xy, x y , y x , 1 xy respectively, see Figure 4. As both objective function and constraints are symmetric with respect to permuting x and y, the minimum of f {y} in O y is same as the minimum of f {x} in O x . Therefore, it suffices to minimize over three regions instead of all four.

For Problem (2) specifically, this means: Proposition 3.2 There exists a natural action of the symmetric group S 4 on the space of molecules (R >0 ) 2 4 by permuting the 4 binding sites:

S 4 × (R >0 ) 2 4 -→ (R >0 ) 2 4 , (σ, (w I )) -→ σ • (w I ) := (w σ(I) ).
This action permutes the regions O I and preserves binding polynomials as well as absolute interactions, i.e., Φ(σ

•W) = Φ(W) and σ •W = W for all W ∈ (R >0 ) 2 4 and all σ ∈ S 4 .
Instead of computing in all 2 11 regions, it therefore suffices to consider one region per S 4 -orbit of regions. However, this still requires 249 regions to be checked, which is still too much for practical purposes.

3.3.

Reducing the degree of the constraints. Given a binding polynomial P , the polynomial equations in Problem (3) arising from its coefficients a 1 , . . . , a 4 pose a serious computational problem. To sidestep this issue we now introduce new coordinates s I := I ⊆I w I . Here, s I represents the free energy difference of microstate I to the fully unoccupied state. On the positive orthant (R >0 ) 2 n , this defines a bijection

(R >0 ) 2 n (R >0 ) 2 n (w I ) W = I ⊆I w I (s I ) = S I ⊆I s (-1) |I\I | I ϕ ϕ -1 ∈ ∈ ∈ ∈
In the new coordinates, the formerly polynomial constraints a k = |I|=k I ⊆I w I become linear constraints a k = |I|=k s I for k = 1, . . . , n. While the functions f I become more complicated in the new coordinates, they still remain monomial and thus it has minimal impact on the runtime: For example, for I = {{1, 2, 3}}

f I (W ) = w -1 123 • |I|>1, I ={1,2,3} w I becomes f I • ϕ -1 (S) = s 1234 s 2 12 s 2 13 s 2 23 s 2 123 s 3 1 s 3 2 s 3 3 s 4 .
In the new coordinates, Problem (3) becomes: minimize (5)

(w I ) I⊆[4] r subject to f I • ϕ -1 (S) ≤ r for all I ⊆ {I ⊆ [n] | |I| > 1}, a 1 = s 1 + s 2 + s 3 +
3.4. Solving the polynomial optimization problem. We will use SCIP [START_REF] Gleixner | The SCIP Optimization Suite 6.0[END_REF] to solve polynomial optimization Problem (5), which is currently one of the fastest non-commercial solvers for non-linear programming. It uses a branch and bound method to solve the optimization problem with non-linear constraints. Example 3.3 shows a computation using SCIP for a molecule with 3 binding sites.

Example 3.3 Consider the polynomial P 1 of Example 2.7. Figure 5 shows the full input on the left and the partial output on the right. In the input, the first constraints c1, c2, c3 enforce |I|=k s I = a k for i = 1, 2, 3. The remaining constraints c4 to c19 enforce f I • ϕ -1 (S) ≤ r. For example, c19 states that s 1 s 2 s 3 /s 123 ≤ r which is equivalent to f ∅ (W ) = (w 123 w 12 w 13 w 23 ) -1 ≤ r in the coordinates w I . The output states an approximate optimal value of r = 13.5 and lists all values of s I of the minimal molecule, which gave rise to the numbers in Example 2.7. Unfortunately, the 2 11 constraints in Problem (5) make it practically impossible to solve directly. Instead, we focus on computing an upper bound b + and a lower bound b -for the minimal interaction. Both bounds are easier to compute, because they arise from problems with significantly less constraints. Our aim is to obtain bounds which either coincide or are sufficiently close to each other to allow us to compare the different minimal interactions with each other.

0 < s1 <2 0 < s2 <2 0 < s3 <2 0 < s12 <3 0 < s13 <3 0 < s23 <3 4 <= s123 <=4 1 <=
The upper bound b + . Our upper bound b + for the solution to Problem (5) is the solution to Problem (4) for a region O I . In other words, instead of minimizing the interaction over all molecules in (R >0 ) 2 4 , we only minimize it over the molecules in O I . It is an upper bound as the global minimum need not be in the region that we have chosen. Computing the solution to Problem (4) still required changing to the new coordinates (s I ). Our region of choice was O ∅ , i.e., the region where w I ≥ 1 for all I ∈ [4] >1 , as it was the only region in which all computations terminated successfully during testing.

The lower bound b -. Our lower bound b -is the solution to a relaxed Problem (5) with some of the 2 11 constraints of the form f I • ϕ -1 (S) ≤ r removed. It is a lower bound since removing constraints enlarges the feasible set over which the minimum is computed. Appendix A contains a list of constraints that remain in the final computation. We picked a region O I in each S 4 orbit of regions and only considered constraints given by functions f I that are dominant in one of the regions.

Experimental Results

In this section, we present the results of computational experiments on data from [START_REF] Imai | Analyses of oxygen equilibriums of native and chemically modified human adult hemoglobins on the basis of Adir's stepwise oxygenation theory and the allosteric model of Monod, Wyman, and Changeux[END_REF][START_REF] Ikeda-Saito | Thermodynamic properties of oxygen equilibria of dimeric and tetrameric hemoglobins from Scapharca inaequivalvis[END_REF] which is also summarized in [START_REF] Connelly | Analysis of zeros of binding polynomials for tetrameric hemoglobins[END_REF]. Section 4.1 summarizes the origins of the data, Section 4.2 contains the bounds on the minimal interactions, and Section 4.3 contains molecules attaining the upper bound. 4.1. Data. The first data set consists of eight oxygen binding polynomials (P 1 to P 8 ) of human adult hemoglobin in its native form and chemically treated with iodoacetamide, N-ethylmaeimide, and carboxypeptidase A. The oxygen binding was observed both in absence and in presence of 2,3-diphosphoglycerate (DPG) [START_REF] Imai | Analyses of oxygen equilibriums of native and chemically modified human adult hemoglobins on the basis of Adir's stepwise oxygenation theory and the allosteric model of Monod, Wyman, and Changeux[END_REF]. One of the goals was to determine how the chemical treatments and DPG affect the binding affinity of hemoglobin, the results are indicated by in Table 1.

The second data set contains five binding polynomials (P 9 to P 13 ) of native hemoglobin HbII of Scapharca inaequivalvis measured at different temperatures [Ike+83, Table 3], see Table 2. One of the goals was to investigate how temperature affects cooperativity in oxygen binding.

Table 3 lists the coefficients a 1 , a 2 and a 3 of the binding polynomials P 1 to P 13 , the coefficients a 0 and a 4 are normalized to 1. Furthermore, it lists the maximal slope to the Hill plot n max , as reported by [START_REF] Imai | Analyses of oxygen equilibriums of native and chemically modified human adult hemoglobins on the basis of Adir's stepwise oxygenation theory and the allosteric model of Monod, Wyman, and Changeux[END_REF] and [START_REF] Ikeda-Saito | Thermodynamic properties of oxygen equilibria of dimeric and tetrameric hemoglobins from Scapharca inaequivalvis[END_REF], which relates the variance of the probability distribution on the macrostates {0, 1, 2, 3, 4} at the respective ligand activity to the variance of a binomial distribution with the same mean [START_REF] Abeliovich | An empirical extremum principle for the hill coefficient in ligand-protein interactions showing negative cooperativity[END_REF][START_REF] Martini | Cooperative binding: a multiple personality[END_REF]. A value larger than 1 indicates a variance higher than the maximal value an independent system can generate. Table 2. Temperature in degree Celsius at which binding polynomials 9 to 13 were determined.

Minimal absolute interactions.

As explained in Section 3.4, solving Problem (5) for n = 4 is very hard, which is why we computed bounds instead. The upper bound is the minimal absolute interaction in the region O ∅ , i.e., the minimal interaction in the region where all interactions are 1 or larger. The lower bound comes from a relaxation of Problem (5) which enlarges the set of molecules over which the minimum is computed. Looking at the bounds b + and b -in Table 3, we observe two things:

(1) Apart from P 2 , P 3 , and P 8 , the bounds determine the minimal interaction up to a precision of 15%. (2) The bounds are almost sufficient to rank the minimal interactions of P 1 to P 8 and P 9 to P 13 : P 6 > P 4 > P 1 > P 3 > P 5 > P 7 > P 8

and P 9 = P 10 = P 12 > P 11 > P 13

with P 2 possibly being larger than P 6 or between P 4 and P 1 .

Comparing the ranking of the binding polynomials by their minimal interactions with the ranking by the maximal slope of their Hill plot n max reveals the following similarities and differences:

P 9 to P 13 . We begin by looking at polynomials P 9 to P 13 as there are prominent differences between minimal interactions and maximal slopes of the Hill plots: While all maximal slopes are approximately equal, there are stark differences in the minimal interactions. This is because the minimal interaction depends on the entire binding curve, while the maximal slope does not. Figure 6 illustrates that while the binding curves of P 9 to P 13 may have identical maximal slopes in logarithmic coordinates, there are subtle differences between both curves. The binding curve of P 11 has a stronger sigmoidal shape, which is why more interaction is necessary to generate the binding behaviour. P 1 to P 8 . While the ranking by minimal interaction is different to the ranking by maximal slope, we retain nearly all important relations in Table 1, which indicate how the chemicals affect the binding behaviour. The only exceptions are P 4 and P 6 , where P 4 has a higher maximal slope and P 6 has a higher minimal interaction. This is not surprising, as Figure 7 illustrates that the binding curve P 4 has a higher maximal slope, while the binding curve P 6 has a stronger curvature. Another explanation can be found in the coefficients of P 4 and P 6 : In P 6 , coefficient a 2 is larger than a 1 and a 3 and almost same size as a 4 = 1. This means that the distribution on the macrostates stabilizes around macrostate 2. This binding behaviour does not only requires more interaction between the binding sites, it also prevents large variances of the distribution on the macrostates and thus leads to a reduction of n max . In P 4 , coefficients a 1 , a 2 , a 3 are smaller than a 4 = 1. This means that the distribution on the macrostates leans toward macrostate 4, which leads to a higher value of n max .

a 1 a 2 a 3 b + b -n max P 1 0.

Molecules solving the optimization problem. Table 4 lists all molecules

realizing the upper bound b + of the minimal interaction. We see that the many interaction energies are close to 1. This is not surprising since we are minimizing the exponential of an L 1 -norm and L 1 optimizations are known to provide sparse solutions [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

Looking at the structure of the solutions in more detail, we see that most solutions include a non-trivial weight for a pairwise interaction (w 34 ). The only polynomial whose molecule possesses trivial pairwise interaction is P 2 which has the highest n max . This is plausible as a high n max means a high variance of the macroscopic distribution, which in turn implies more weights towards macrostates 0 and 4. And any nontrivial pairwise interaction would put more weight on macrostate 2.

Moreover, we mentioned how the macroscopic distribution for P 6 stabilizes around 2 while the macroscopic distribution for P 4 leans towards 4. This can be seen in their minimal molecules in O ∅ , as the molecule for P 6 has large w 34 and the molecule for P 4 has a large w 1234 .

Discussion and Outlook

5.1. Cooperativity and minimal absolute interaction. While cooperativity has always been a property of the macroscopic binding behavior, it is often interpreted on a microscopic level, for instance that "binding of a ligand molecule increases the receptor's apparent affinity, and hence increases the chance of another ligand molecule binding" [START_REF] Stefan | Cooperative binding[END_REF]. There is however no straightforward connection between the macroscopic and microscopic properties, as generally there are infinitely many different molecules with the same binding curve. Thus creating a bridge between macroscopic and microscopic properties usually relies on additional restrictions, such as having indistinguishable sites.

In this paper, we considered the minimal interaction necessary to generate a fixed overall binding curve. In other words, we study the minimal deviation from a molecule with no interaction that is needed to create a prescribed binding behavior. We showed that computing the minimal interaction gives rise to a challenging polynomial optimization problem, and proposed easier ways for obtaining upper and lower bounds. Applying our techniques to overall binding curves and comparing our results to the maximal slope of the Hill plot, we observe:

(1) The upper and lower bounds were sufficient to rank all but one of the curves.

(2) The complete ranking of the curves by minimal interaction may differ from the ranking by maximal Hill slope, but most of the important relations are the same.

(3) In the cases were minimal interaction and maximal Hill slope disagreed, the curve with the higher minimal interaction displayed a stronger curvature (see Figures 6 and7).

Conceptually, the minimal interaction considers the entire curve and captures any deviation from an independent system, while the maximal Hill slope only indicates an abnormally high variance of the macroscopic distributions. This is also reflected in our computations, as our results based on the minimal interaction largely coincides with the results based on the maximal Hill slope and the small differences in the results could be explained by signs of interaction outside the region of maximal Hill slope. Additionally, the corresponding minimal molecules may give valuable insights into the microscopic dependencies of the molecule. While the minimal molecule may not be a physically correct representation of the actual chemical molecule that was measured, in a way it describes the simplest possible set of interactions required for the binding curve. Note that we chose to minimize the L 1 -norm of the actual interaction energies because it was the most natural norm in our framework. Alternatively, other L p -norms could be used, and it is unclear whether they yield better qualitative results. However, with current methods the minimal interaction is much harder to compute than the maximal Hill slope and thus new mathematical approaches are required. 5.2. The mathematical problem. Ligand binding raises many intriguing questions in applied algebra and geometry, and stands to benefit equally from further insights into the intrinsic geometry of the problems. While previous works have explored the application of polynomial system solving techniques to find decoupled molecules in a situation with different types of ligands [RMT19; MSU13], this work explored the use of algebraic optimization to calculate the minimal interaction required to generate a certain binding behavior.

We illustrated how searching for the minimal interaction energy required to generate a binding curve leads to a problem in algebraic optimization. Even though, computational tools are readily available and the structure of the problem seems relatively simple on first sight, it becomes complicated when more than 3 binding sites are considered. In particular the fact that we had to calculate upper and lower bounds by a relaxed problem and were not able to close their gap completely, shows that our approach is not yet ready for a simple ready-to-use tool to quantify cooperativity. Even though we were able to rank the polynomials in Section 4.2, we were not able to compute precise minimal energies for all polynomials. We thus believe further research into the computational problem is warranted, though it may need to come from a different direction.

Exploiting symmetries. One possibility would be to show that the minimal molecule have some symmetry, which would reduce the number of variables (interaction energies) in our problem. While both Problem 5 and the molecules in Table 4 show some kind of symmetry, it is mathematically not clear whether that always needs to be the case. There are results using representation theory which relate symmetries in polynomial optimization problems to symmetries in their solutions [START_REF] Riener | Real root finding for equivariant semialgebraic systems[END_REF], however the assumptions on the symmetry in the polynomial optimization problem are quite strong and not satisfied by Problem 5.

Tropical Geometry. Tropical geometry [START_REF] Maclagan | Introduction to tropical geometry[END_REF] revolves around functions on the maxplus semiring, such as the absolute interaction in logarithmic coordinates: For our problem, L is the image of the feasible set given by the conditions Φ(W ) = P which is an exponential variety of codimension n. Theoretically, there exist only a finite number of intersection patterns of L and T . Classifying these intersection patterns would allow us to compute minimal interactions and minimal molecules almost instantly.

Tensors. Tensors [START_REF] Landsberg | Tensors: geometry and applications[END_REF] are multidimensional arrays of real numbers which generalize matrices. We can regard molecules as (positive) tensors in R 2ו••×2 =: (R 2 ) ⊗n , whose entries are indexed by {0, 1} n which is in bijection with the subsets of [n].

From this point of view, many concepts in tensor theory have natural analogues in ligand binding, such as symmetric tensors and molecules with indistinguishable sites. Notably, a rank one tensors is up to scaling a tensor with entries (1, 0, . . . ) : s 1 , (0, 1, 0, . . . ) : s 2 , (1, 1, 0, . . . ) : s 1 s 2 , . . . for some s i ∈ R.

Hence rank one tensors correspond to molecules with independent sites. Therefore, decomposing a given tensor into rank one tensors translates to decomposing a given molecule into different molecules with independent sites whose mix has the same binding curve. And the tensor rank translates to the minimal number of different molecules with independent sites needed to recreate the same binding curve. Thus the tensor rank of a molecule could serve as an interesting, discrete measure for cooperativity.

Dependency measures. From a stochastic or statistical point of view, the problem of cooperativity is a problem of quantifying (minimal) dependency, which is a classical research topic in these fields [Rén59; SW81; Koy+87]. In particular copulas [Sch91; DS10; Nel07] have become a common tool to describe dependencies between random variables. Copulas cannot be applied directly to our problem since they are defined on continuous, not discrete, random variables. Moreover, they aim at modeling the dependency of given random variables which in our scenario corresponds to the dependency of the sites of a known molecule. A copula then relates the one-dimensional marginal distributions to the multivariate minimization problem of finding the minimal dependency required to generate a given distribution of a sum of dependent variables is usually not treated. Nevertheless, the problem of cooperativity should also be treated from a more probabilistic point of view and the available concepts should be considered in more detail in the context of ligand binding.

Figure 1 .

 1 Figure 1. The oxygen concentration in dog and horse blood in relation to the ambient oxygen concentration [BHK04, Fig. 2].

interaction energy w 12 w 134 Figure 2 .

 1342 Figure 2. A molecule with 4 binding sites.

  also commonly known as the Adair equation [ABJ25; SL13], with positive real coefficients a k := |I|=k I ⊆I w I ∈ R >0 .

Figure 3 .

 3 Figure 3. Minimal molecules for binding polynomial P 1 , P 2 , P 3 .

Figure 4 .

 4 Figure 4. Partitions in plane

  -r * s1 * s2 * s3 <=0 c5 : s1 * s2 * s123 -r * s3 * s12 ^2 <=0 c6 : s1 * s3 * s123 -r * s2 * s13 ^2 <=0 c7 : s2 * s3 * s123 -r * s1 * s23 ^2 <=0 c8 : s12 ^2* s13 ^2* s23 ^2 -r * s1 ^3* s2 ^3* s3 ^3* s123 <=0 c9 : s1 ^3* s2 * s3 * s123 -r * s12 ^2* s13 ^2 <=0 c10 : s1 * s2 ^3* s3 * s123 -r * s12 ^2* s23 ^2 <=0 c11 : s13 ^2* s23 ^2 -r * s1 * s2 * s3 ^3* s123 <=0 c12 : s1 * s2 * s3 ^3* s123 -r * s13 ^2* s23 ^2 <=0 c13 : s12 ^2* s23 ^2 -r * s1 * s2 ^3* s3 * s123 <=0 c14 : s12 ^2* s13 ^2 -r * s1 ^3* s2 * s3 * s123 <=0 c15 : s1 ^3* s2 ^3* s3 ^3* s123 -r * s12 ^2* s13 ^2* s23 ^2 <=0 c16 : s1 * s23 ^2 -r * s2 * s3 * s123 <=0 c17 : s2 * s13 ^2 -r * s1 * s3 * s123 <=0 c18 : s3 * s12 ^2 -r * s1 * s2 * s123 <=0 c19 : s1 * s2 * s3 -r * s123 <=0 Bounds

Figure 5 .

 5 Figure5. Computing the minimal absolute interaction for a molecule with 3 sites using SCIP.

Figure 6 .

 6 Figure 6. The binding curves of P 9 (black) and P 11 (red) in normal (left) and logarithmic (right) scale.

Figure 7 .

 7 Figure 7. The binding curves of P 4 (black) and P 6 (red) on [0, 4] (left) and zoomed in on [0, 0.2] (right). Note how the red curve starts and ends lower as the black curve, but overtakes it in between.

(

  R >0 ) 2 n R 2 n (log(w I )) =: (x I ) (w I ) I⊆[n] max(w I , w -1 I ) I⊆[n] max(x I , -x I ) =: (x I ) trop ∈ ∈In tropical geometry, the function • trop is called a tropical rational function. It is a piecewise-linear function, linear on the image of each O I , and the boundary separating its regions of linearity is called its tropical hypersurface T . Since the gradient of • trop is constant on each region of linearity, optimizing • trop on any set L boils down to understanding the intersection of L with T .

  w 1 w 2 w 12 + w 1 w 3 w 13 + w 1 w 4 w 14 +w 2 w 3 w 23 + w 2 w 4 w 24 + w 3 w 4 w 34 , a 3 = w 1 w 2 w 3 w 12 w 13 w 23 w 123 + w 1 w 2 w 4 w 12 w 14 w 24 w 124 +w 1 w 3 w 4 w 13 w 14 w 34 w 134 + w 2 w 3 w 4 w 23 w 24 w 34 w 234 , a 4 = w 1 w 2 w 3 w 4 w 12 w 13 w 14 w 23 w 24 w 34 w 123 w 124 w 134 w 234 w 1234 .

  w 1 w 2 w 12 + w 1 w 3 w 13 + w 1 w 4 w 14 +w 2 w 3 w 23 + w 2 w 4 w 24 + w 3 w 4 w 34 , a 3 = w 1 w 2 w 3 w 12 w 13 w 23 w 123 + w 1 w 2 w 4 w 12 w 14 w 24 w 124 +w 1 w 3 w 4 w 13 w 14 w 34 w 134 + w 2 w 3 w 4 w 23 w 24 w 34 w 234 , a 4 = w 1 w 2 w 3 w 4 w 12 w 13 w 14 w 23 w 24 w 34 w 123 w 124 w 134 w 234 w 1234 .

Table 1 .

 1 Table 3 lists the upper (b + ) and lower (b -) bounds for the minimal absolute interaction for the binding polynomials P 1 to P 13 . Conditions under which binding polynomials 1 to 8 were measured. P i ≺ P j means P i that the Hill plot of P i has a lower maximal slope compared to P j . Clam HbII at 10 • at 15 • at 20 • at 25 • at 30 •

	Human Hb	in the absence of 2mM DPG	in the presence of 2mM DPG
	untreated	P 1	≺		P 2
	treated with iodoacetamide	P 3	≺		P 4
	treated with N-ethylmaeimide	P 5	≺		P 6
	treated with carboxypeptidase A	P 7	≺		P 8
	P 9	P 10	P 11	P 12	P 13

Table 3 .

 3 Binding polynomials, bounds on the minimal absolute interaction and n max as reported by[START_REF] Ikeda-Saito | Thermodynamic properties of oxygen equilibria of dimeric and tetrameric hemoglobins from Scapharca inaequivalvis[END_REF][START_REF] Imai | Analyses of oxygen equilibriums of native and chemically modified human adult hemoglobins on the basis of Adir's stepwise oxygenation theory and the allosteric model of Monod, Wyman, and Changeux[END_REF]. Coefficients a 0 and a 4 are normalized to 1 for all polynomials.

		835 0.379 0.541 527 527 2.51
	P 2 0.789 0.154 0.0648 3322 1600 3.09
	P 3	1.42 2.42 0.752 111	63 1.63
	P 4 0.647 0.568 0.0986 2002 1460 2.71
	P 5	2.0 2.31	2.04	16	16 1.44
	P 6 0.539 0.909 0.554 3033 3030 2.27
	P 7	3.47 4.74	2.76 2.27 1.68 1.15
	P 8	3.26 5.36	2.23 7.64 2.19 1.23
	P 9	1.4	1.0	0.62	66	66 2.08
	P 10	1.4 0.96	0.60	66	66 2.10
	P 11	1.2 0.93	0.70 123 123 2.08
	P 12	1.4 0.95	0.62	66	66 2.10
	P 13	1.1 0.98	0.59 175 175 2.12

Table 4 .

 4 12 • • • w 24 w 34 w 123 w 124 w 134 w 234 w 1234 Molecules of P 1 to P 13 realizing the upper bound b + . All interactions strictly larger than 1 are highlighted in blue.

	w 1 w 4 w P 1 0.209 0.209 0.209 0.209 w 2 w 3	1.00	3.70 1.00 1.00 14.5 1.00 9.79
	P 2 0.0822 0.0822 0.0822 0.543	1.00	1.00 6.60 1.00 14.7 1.00 34.3
	P 3 0.179 0.179 0.531 0.531	1.00	7.12 1.00 1.00 1.00 1.00 15.6
	P 4 0.101 0.101 0.223 0.223	1.00	9.42 1.00 1.00 1.00 1.00 212
	P 5 0.500 0.500 0.500 0.500	1.00	4.24 1.00 1.00 1.69 1.69 1.32
	P 6 0.135 0.135 0.135 0.135	1.00	45.1 1.00 1.00 2.49 2.49 10.9
	P 7 0.569 0.569 1.17 1.17	1.00	1.30 1.00 1.00 1.00 1.00 1.75
	P 8 0.265 0.265 1.36 1.36	1.00	2.06 1.00 1.00 1.00 1.00 3.70
	P 9 0.350 0.350 0.350 0.350	1.00	3.16 1.00 1.00 1.97 1.97 5.43
	P 10 0.350 0.350 0.350 0.350	1.00	2.84 1.00 1.00 2.11 2.11 5.26
	P 11 0.300 0.300 0.300 0.300	1.00	5.33 1.00 1.00 2.24 2.24 4.60
	P 12 0.350 0.350 0.350 0.300	1.00	2.76 1.00 1.00 2.26 2.26 4.73
	P 13 0.275 0.275 0.275 0.275	1.00	7.96 1.00 1.00 1.66 1.66 8.01

Acknowledgments. The authors would like to thank Felipe Serrano for valuable technical support with SCIP. The first and second author were partially supported by the Max Planck Institute for Mathematics in the Sciences.

Yue Ren is supported by the UKRI Future Leaders Fellowship (MR/S034463/1). Mohab Safey El Din is supported by the PGMO grant CAMiSAdo and the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement N. 813211 (POEMA).

Email address: jwrmartini@gmail.com


Samuel I Billong

G E Kouamou

T Bouetou
email: tbouetou@gmail.com

Bruce Watson

Eric Badouel

Eds. Oumar Niang

Discretization of an adaptive prediction model using a non-homogeneous meshing

Keywords: Prediction model, Diffusion equation, Complex geometry, discretization, finite volume methods

In a previous work, we have developed a model based on 2D diffusion equations which was discretized on homogeneous meshes taking into account the mobility of entities in an existing prediction model. This model does not fit with the reality of administrative boundaries. The purpose of this paper is to deal with a more realistic approach where the boundary of the study regions are irregular. The work is done in a meshing imposed by the geographical configuration of the administrative divisions and the discretization is done with the finite volumes method for complex geometries. The validation is done using a set of administrative divisions through a multi-agent system and a shapefile extracted in a geographic information system were the modified SIR epidemiological model is applied to highlight the dynamic effect in the spread.

Introduction

The numerical evaluation and implementation of a continuous model is possible through the discretization, which is the first step in the process of transferring the equations into the discrete counterparts [START_REF] Ciarlet | [END_REF][START_REF] Vignal | Numerical and theoretical study of a dual mesh method using finite volume schemes for two phase flow problems in porous media[END_REF]. This transformation usually simplify some aspects of the reality as well as it approximates the hypothesis of the models. In the prediction domain like other fields of research, they are some cases where the developed theories do not confirm the reality. In this work we have gone through several prediction model in epidemiology [3, 4, 8, and 9]. The gap between the simulation and the reality of the field is due to the assumption that are made during the transformation of the model. Concerning the domain of complex systems, particularly the prediction of phenomena, the dynamic aspect of mobile entities which are moving in the studied environment is not consider.

In the previous studies [START_REF] Billong | A spatio-temporal Model for phenomena dynamics based on 2D Diffusion equations[END_REF], a diffusion equation was written to manage the mobility of the entities. Since their displacement is random, a probabilistic matrix was introduced to integrate the adaptation of the model. In the construction of this model, which resulted in 2D diffusion equations, approximations have been made. With all these assumptions, it is clear that the model implemented does not represent the reality reliably then it would be normal that certain results do not comply with the observation as we specified above. Based on these remarks, we come up with the idea of making the model more realistic by first changing the forms of the sites to get closer to the real geographical forms.

The purpose of this article is to improve the discretization of the established model in order to fit closer to reality through the application of this model independently to the geometries of the studied region. To this effect, the initial discrete model obtained with the finite volumes method on regular meshes is modify considering non homogeneous meshes to manage complex geometries. This approach is due to the real configuration of the study sites whose shapes are mainly arbitrary. For illustration, four sites in a real environment are chosen from the administrative divisions through a GIS to constitute a closed environment where simulations will be made.

The rest of this paper is structured as follows. The section 2 presents the literature review where the basic theories and the former development are described. The section 3 presents the methodology, the manner in which the model is built with finite volume discretization for complex geometries. Section 4 describes a hybrid SIR obtained by applying this idea on the traditional SIR model. We will end in section 5 with a conclusion and further works.

Background.

Population dynamics and particles diffusion

The foundation of this study is focused on the equations taken in [START_REF] Garnier | Analyse mathématique de modèles de dynamique des populations : équations aux dérivées partielles paraboliques et équations intégro-différentielles[END_REF] with the following forms

𝜕𝑢(𝑡,𝑥) 𝜕𝑡 = 𝐷(𝑢)(𝑡, 𝑥) + 𝑓(𝑥, 𝑢(𝑡, 𝑥)) 𝑡 > 0 , 𝑥 ∈ 𝑅 (1)
These equations occur in a wide variety of fields such as combustion, chemistry, biology or ecology. They are used to model the evolution of the entities that interact with each other and move. Applied to the population dynamics, the quantity u (t, x) represents the population density at the period t at the position x. The reaction term f (x, u) corresponds to the growth rate of the population. The movement of individuals is described by the dispersion operator D. Depending on the mode of movement of the individuals, this operator will be local or non-local. Our attention is focused mainly on a single type of reaction-dispersion equation where the dispersion operator D is an elliptic differential operator of the second order.

𝜕𝑢(𝑡,𝑥) 𝜕𝑡 = 𝜕 2 𝑢(𝑡,𝑥) 𝜕𝑥 2 + 𝑓(𝑥, 𝑢(𝑡, 𝑥)) (2)
In this first family of equations we found that the equation of dimension is not valid, equation is purely mathematical but the objective of this work is to deal with the quantities reliable with the reality. This motivated a balance approach that originated in the particles diffusion In [START_REF] Billong | A spatio-temporal Model for phenomena dynamics based on 2D Diffusion equations[END_REF] it is found that the equations resulting from the particles diffusion come from the balances of mobile entities in a control volume, this approach was built with the help of [START_REF] Salamito | Damien physique tout en un Editeur[END_REF].

In the same way, [START_REF] Wu | Spatial Structure: partial differential equations models[END_REF] starting with the basic concepts developed a model for the spatial spread of diseases involving hosts in random displacement during certain stages of the progression of the disease. It led to a diffusion model based on the conservation law and Fick's law. He applied this model to the study of two cases, namely the spread of rabies in continental Europe during the period 1945-1985 and the rate of spread of West Nile virus in North America.

The former development

With the basic equation

𝜕𝑢(𝑡,𝑥) 𝜕𝑡 = 𝐷 ( 𝜕 2 𝑢(𝑡,𝑥) 𝜕𝑥 2 + 𝜕 2 𝑢(𝑡,𝑦) 𝜕𝑦 2 ) (3)
We considered a closed environment with homogeneous site distribution, then the discrete model was suitable. The choice fell on the finite volume discretization method which is adapted to the equations that represent conservation laws [START_REF] Guinot | Méthodes numériques Appliquées (résolution numérique des équations différentielles de l'ingénieur) Polytech[END_REF], [START_REF] Guichard | Schémas volumes finis sur maillages généraux en milieux hétérogènes anisotropes pour les écoulements polyphasiques en milieux poreux Thèse de doctorat Université de Paris -Est novembre[END_REF], [START_REF] Martin | Elaboration de solveurs Volumes Finis 2D/3D pour résoudre le problème de l'élasticité linéaire[END_REF] and [START_REF] Silva | Méthodes et Analyse Numériques[END_REF].

Then the D coefficient were built according to hypotheses. The coefficient D is a probabilistic contribution matrix. The general term of D is 𝐷 𝑖,𝑗 = 𝜔𝑃 𝑖,𝑗 where 𝜔 ≈ (∆𝑥^2)/∆𝑡 depends on the average speed characteristic of the flow moving of mobile entities and 𝑃 𝑖,𝑗 the probabilistic matrix of movements. After applying the finite volume method and assumptions the following numerical scheme were obtained

Methodology of a new approach in real context

The sample regions are approximated by the polygons as shown in Figure 1. 𝑢 𝑛+1 represent the number of mobile entities in the site at the time t + dt and 𝑢 𝑛 is the value at the period time t.

The finite volume method for complex geometry applied to the basic equation using STOKES' formula [START_REF] Silva | Méthodes et Analyse Numériques[END_REF]  ∆𝑡 is the time step in the numerical scheme

The boundary conditions are implicitly considered because the quantities above ( 𝑢 𝑛 and 𝑢 𝑛+1 ) are zero everywhere else except on the boundaries between our experimental environments.

Validation and interpretation

The case study is taken from the epidemiology considering the traditional SIR (Susceptible, Infected, and Retired) model [START_REF] Brauer | Mathématical Epidemiology, eds mathematical biosciences subseries[END_REF]. The diffusion is considered through propagation dynamic of the infection.

Figure 2. Transition state diagram of an individual

Note on the Figure 2 that an individual before being in a Remove state (R) must previously passed through the infected state; and only susceptible individuals can become infected. Three administratives regions of the Cameroon are concerned in the illustration denoted H, N, P, R as shown in Figure 3.

Figure 3. Polygonal approximation of the selected regions

In the case of the dynamics of the phenomenon taking into account the spatiotemporal dynamics, we have the following discrete equations for the site H knowing that the value u(t, x) is considered for each state S, I, R. These equations are replicated for each regions considering its neighborhood.

{ 𝑆

Simulation

The software tool is built on REPAST SYMPHONY which is a multi-agent development environment. The regions are extracted from a GIS, then the corresponding shapefiles are integrated in REPAST. The equations are encoded as the behavior of the agents. The parameters fixed before running the simulation are mentioned in Table 1. Initial parameters of the simulation Table 1. The probability of leaving an area 0.5

The probability of changing the state from S to I 0.5

The probability of changing the state from R to S

0.05

The diffusion velocity 10km/days

Discussion

Starting from the differential equations of the propagation of an epidemic in SIR model, we obtained a hybrid equations system coupled to the model developed to manage the dynamics of the mobile entities. By analyzing the comparative histogram (see the annexes), the following observations are made:

 High density areas carry the epidemic faster. Having the same number of agents, on all the sites, it is the smallest sites in area which have a high rate of contamination, here: H and K;  The epidemic persists longer when there is diffusion. This is the case of the region N where the density is low. The epidemic, when there is no diffusion disappears after a few days for lack of contact. Yet the diffusion shows that patients from neighboring areas are reviving the epidemic after the infected natives in the area are cured.  It can be concluded from the analysis of the data presented that the diffusion of mobile entities has an impact on the life cycle of an epidemiological phenomenon also in the case of complex geometry. This observation validates the hypothesis make at the beginning of this work.

Conclusion

The study carried in this paper focus on the discretization on non-homogeneous meshes of an adaptive prediction model built on the basis of diffusion equations. The approach adopted is based on the finite volume method for complex geometries for a mesh imposed by the geometry of the study environment which can be mapped to the administrative division of a territory.

First, it was necessary to recall the basis of the model. Then, the hypothesis were extended to meet the resulting equations from the discretization in a non-homogeneous mesh implemented in the agents as their behavior for simulation. For cartographic visualization, we extracted the shapefile from the sample regions using a GIS. We integrated this shapefile into the multi-agent system were an implementation of the resulting model has been made using the software REPAST SYMPHONY. The observation shows that the consideration of spatiotemporal dynamics in an SIR epidemiological model presents a better reality than the traditional SIR.

The future works will focus on the experimentation in a live situation in order to release the accuracy of the substantial contribution, the improvement of the discretization on any meshing which does not relate to the center of gravity.

  𝑢 𝑖+1,𝑗,𝑥 𝑛 𝛿 𝑙,𝑥 + 𝑢 𝑖-1,𝑗,𝑦 𝑛 𝛿 𝑙,𝑦 ) + ∆𝑡𝐷 𝑙,𝑘 ∆𝑦 2 (𝑢 𝑖,𝑗+1,𝑧 𝑛 𝛿 𝑙,𝑧 + 𝑢 𝑖,𝑗-1,𝑡 𝑛 𝛿 𝑙,𝑡 ) + (1 )𝑢 𝑖,𝑗,𝑘 𝑛 𝛿 𝑙,𝑘

Figure 1 .

 1 Figure 1. The sample regions and the polygonal approximations The following considerations are made:  XH and XN are the position of center of gravity of two any regions denoted by (H) and (N)  |𝜎 𝐻𝑁 | the length of the ridge that separates H and N perpendicularly to the line (XHXN)  |𝐻| measurement of control volume of the region H  𝑑 𝐻𝑁 measurement of the distance XH -XN  𝜕𝐻 set of all the outline of any region H Let remember the basic equation of the model (4)

Figure 4 .

 4 Figure 4. Infection at the initialization, then after 5, 10 and 20 cycles respectively

  in any control volume H, leads to

	𝑢 𝐻 𝑛+1 = ∑	𝑁 ∈ 𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑢𝑟(𝐻)	𝐷 𝑁𝐻 |𝜎 𝐻𝑁 |∆𝑡	𝑢 𝑁 𝑛 -𝑢 𝐻 𝑛 |𝐻|𝑑 𝐻𝑁	+ 𝑢 𝐻 𝑛	(6)
	 𝑢 𝐻

𝑛+1

represents the number of mobile entities at the time t + dt in H  𝐷 𝑖,𝑗 = 𝜔𝑃 𝑖,𝑗 represent an element of the diffusion matrix where 𝜔 is the average speed of diffusion and 𝑃 𝑖,𝑗 the matrix of random contributions between sites.

Table 1 . Initial parameters of the simulation

 1

	Parameter	Value	Parameters		value
	Susceptible agent per region	200	The	distance	of	500
			contamination		meters
	Infected agent per region	1	The	probability	of	0.2
			changing the state from I to
			R
	The number of agents retired	0
	per zone

Proceedings of CARI 2020

Discretization of an adaptative prediction model using a non-homogeneous meshing