Fabrizio Clarelli
email: clarelli@iac.rm.cnr.it

Cristiana Di Russo
email: dirusso@math.univ-lyon1.fr

Roberto Natalini
email: roberto.natalini@cnr.it

Magali Ribot
email: ribot@unice.fr

A fluid

INTRODUCTION

A biofilm is a complex gel-like aggregation of microorganisms like bacteria, cyanobacteria, algae, protozoa and fungi. They stick together, they attach to a surface and they embed themselves in a self-produced extracellular matrix of polymeric substances, called EPS. Even if a biofilm contains water, it is mainly in a solid phase. Biofilms can develop on surfaces which are in permanent contact with water, i.e. on solid/liquid interfaces or on different types of interfaces such as air/solid, liquid/liquid or air/liquid.

The term "biofilm" was proposed in 1978 by Costerton [START_REF] Costerton | How bacteria stick[END_REF] to describe what is now considered as the preferred form of microbial life in natural, industrial and hospital environments. Some biofilms are useful, providing valuable services to human society or to the functioning of natural ecosystems. Bacteria in the subsurface normally grow as biofilms on the soil matrix, for instance, and can help to remove contaminants from the soil or ground waters. However, other biofilms are harmful, causing serious health and economic problems. Their propensity for attachment causes problems in many inrastructures, like industrial pipelines, ship hulls, nuclear power stations, air conditioning systems, water distribution systems. Hospitals are also susceptible to colonization by microorganisms growing in biofilms [START_REF] Hall-Stoodley | Bacterial Biofilms: from the natural environment to infectious diseases[END_REF]. Harmful biofilms often develop even under adverse conditions so that removing them is often difficult. Since biofilms play a significant role in many natural and engineered systems, understanding the mechanisms of biofilm formation, growth, and removal could be the key for promoting biofilms which have useful applications and also for reducing deleterious biofilms.

Biofilms are not simply passive assemblages of cells, but they form structurally and dynamically complex biological systems. Their development is often characterized as a multistage process. First, some free-floating bacteria approach the surface and get attached. Then, during a phase of colonization, bacteria loose flagella and produce EPS. During the phase of growth, bacteria build a 3D structure, influenced by a variety of environmental factors. In the end, a part of the biofilm may detach itself in order to colonize other parts of the surface.

Since the topic is huge and of great interest, some mathematical models have already been proposed. At the beginning, mathematical modeling of biofilm was mainly focused on predicting growth balance, sometimes with practical applications in mind, as in [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF][START_REF] Wanner | Competition in biofilms[END_REF][START_REF] Wanner | A multispecies biofilm model[END_REF]. These are generally 1-D models with reaction-diffusion equations for nutrients and other substrates, sometimes with a moving boundary. The first multidimensional models were discrete and based on cellular automata. For example, models proposed by the Delft's team [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF] are mainly multidimensional, multispecies and multisubstrates spatially discrete models, which have been solved by individual-based approach or cellular automata. They are quite exhaustive from the biological point of view, at least qualitatively, but not fully satisfactory, because in many cases it is hard to give well based ruled for all the involved individuals, and for the difficulties to simulate large colonies of millions of individuals with the discrete approach. Besides, it is not easy to predict the actual behavior of the solutions. Fully continuous models have also been considered. The recent review by Klapper and Dockery [START_REF] Klapper | Mathematical Description of Microbial Biofilms, to appear[END_REF] focuses mainly on this kind of models, following the idea of treating a biofilm as a viscoelastic material that expands in response to growth-induced pressure. A mathematical justification of these models was proposed by Overgaard [START_REF] Overgaard | Application of variational inequalities to the moving-boundary problem in a fluid model for biofilm growth[END_REF] using variational techniques. Among them, an important class of models, proposed by Alpkvist and Klapper [START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF], is based on a multidimensional and multispecies description, where biofilms are divided into a biomass region and a liquid region separated by an interface computed by moving fronts techniques. However, it is difficult to find a physical law governing the evolution of the interface. More recently, Zang, Cogan, and Wang [START_REF] Zhang | Phase-field models for biofilms II. 2-D numerical simulations of biofilmflow interaction[END_REF] proposed a two phases model, considering the polymer network and the solvent, and analyze numerically the case of detachment under different initial conditions. Another continuous model, including more biological details, was proposed by Anguine, King and Ward [START_REF] Anguige | A multiphase mathematical model of quorum sensing in a maturing Pseudomonas aeruginosa biofilm[END_REF]. It concerns biofilms produced by the Pseudomonas aeruginosa, a bacterium that causes serious infections. In this multispecies PDEs model, four different phases are considered: live cells, dead cells, EPS, and liquid. The influence of nutrients is also taken into account, as well as quorum sensing, one of the various signaling mechanism of cells, and also some different medical treatments, like antibiotics and antiQS drugs. Transport equations are introduced to model the four phases and the behavior of nutrients, antibiotics and antiQS are given by advection-diffusion equations. A common velocity for bacteria, dead bacteria and EPS is assumed, while a different velocity is taken for the liquid. The system is closed with the no-void condition, together with a supplementary relation between the liquid and EPS, namely, a local increment of EPS causes a local increment of liquid. Thanks to these hypotheses, no equation for the velocities is needed, but the model works only in one space dimension.

In this paper, we propose a new approach, based on mixture theory, to the formation and evolution of cyanobacteria biofilms in several space dimensions. We focus on their development on the stone surfaces of ancient monuments, as for instance fountains walls, i.e.: on stone substrates and under a water layer. The colonization of external surfaces of buildings, monuments and archaeological sites by microorganisms creates an unaesthetically appearance of staining of the stone surfaces [START_REF] Crispim | Cyanobacteria and Biodeterioration of cultural Heritage: A Review[END_REF][START_REF] Cuzman | Biodiversity of phototrophic biofilms dwelling on monumental fountains[END_REF]. It produces also extracellular polymeric substances that cause mechanical stresses to the mineral structure inside the pore system. This can lead to the alteration of pore size and distribution, together with changes in moisture circulation patterns and temperature response. Microorganisms may also alter the water permeability of the minerals by the deposition of surfactants. Cyanobacteria, also known as blue-green algae, are a phylum of bacteria that obtain their energy through photosynthesis, thanks to an elaborate and highly organized system of internal membranes. Our model refers essentially to the order of Chroococcales, in particular Gloeocapsa, which is a genus of photoautotrophic bacteria and a prokaryote, which are relevant for stone damage analysis. In this paper, we assume that bacteria remain attached to EPS and therefore that bacteria and EPS have a common velocity. This is the case for Chroococcales cyanobacteria, since they secrete individual gelatinous capsules which totally embed them (see [START_REF] Congestri | Seasonal succession of phototrophic biofilms in an Italian wastewater treatment plant: biovolume, spatial structure and exopolysaccharides[END_REF][START_REF] Pippo | Capsular polysaccharides of cultured phototrophic biofilms[END_REF]). Even if our focus is on phototrophic species, like cyanobacteria, most of the framework we are going to deal with in this paper can be extended to other species and mixed colonies.

Our first goal is to introduce a model which keeps the physical finite speed of propagation of the fronts. Starting from the ideas of the mixture theory [START_REF] Rajagopal | Mechanics of mixtures[END_REF][START_REF] Astanin | Multiphase models of tumor growth[END_REF], we write balance equations which contain the main assumptions coming from biophysical considerations : mass and momentum conservation, influence of light. For simplicity reasons, we will consider here that nutrients and temperature are not limiting factors and that they do not play a significant role; however, it would be simple to add their dependence in further studies. We take particular care to keep the inertial terms in the momentum equations, to guarantee the hyperbolicity of the system and so the finite speed of propagation. Actually, in most of the models coming from the mixture theory approach, as for instance [START_REF] Farina | On Darcy's law for growing porous media[END_REF], these terms are neglected to simplify the analysis and the numerical approximation, therefore producing diffusive terms which stabilize the fluid and prevent possible breakdowns or other instabilities. However, this simplification introduces a non physical infinite speed of propagation in the problem, and makes difficult to study effectively the evolution of interfaces between the solid (biofilm) and the liquid (water) phases. A possible solution is to use moving fronts techniques (see [START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF]), which lead to other analytical and numerical difficulties and a further approximation in the model. Therefore, we keep the inertial terms and we solve the full hyperbolic problem using robust numerical schemes, as the Riemann Solver-free relaxation schemes [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. There are two important differences with respect to a usual hyperbolic system. First, since we are considering a multiphase fluid, it is difficult to deal with regions where one of the phases may vanish. This is usually solved by neglecting these regions, for instance by selecting adapted initial conditions. In a biofilm this choice is not suitable, since it is important to model also the regions made of biofilm alone or liquid alone. It turns out that it is possible to deal with this problem of vanishing phases by using an Implicit-Explicit scheme in the approximation of the source terms (see Section 3 for more details). The second problem arises from the fact that our model is supplemented with a constraint term due to the mass conservation, which implies that the average hydrodynamic velocity of the mixture is divergence free. This constraint is needed to compute the hydrostatic pressure. To enforce the divergence free constraint, we use a fractional step approach similar to the Chorin-Temam projection scheme [START_REF] Chorin | Numerical simulation of the Navier Stokes equations[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF] for the Navier-Stokes equations, with a quite accurate reconstruction of the pressure term.

The plan of the paper is the following. In Section 2, we present a detailed derivation of our fluid dynamics model, with a specific focus in the modeling of the influence of the light on the growth of the biofilms. The numerical approximation scheme is presented in Section 3. Finally, in the last section, we present numerical tests, with the aim of illustrating the power of our approach and also the influence of the various parameters. Our tests are performed in one, two, and three space dimensions, for a domain with no flux conditions. We study the influence of light on the system and how the estimate of the sound velocity of the medium, which however is not really an easy task from the experimental side, affects the final output and in particular the speed of the front.

THE FLUID DYNAMICS MODEL

To describe the complex structure of biofilms, we consider, as in [START_REF] Anguige | A multiphase mathematical model of quorum sensing in a maturing Pseudomonas aeruginosa biofilm[END_REF], four different components: Live cyanobacteria (B), Dead cyanobacteria (D), Extracellular matrix of Polymeric Substances or EPS (E), and Liquid (L). We denote the concentration of biomass by C φ = ρ φ φ , where ρ φ is the mass density of a phase in [g/cm 3 ] and φ = B, D, E, L is the volume fraction of the phases. We assume that the biomasses are incompressible, so that ρ B , ρ D , ρ L and ρ E are positive constants. We also assume that the phases have all the same constant density, which is a good approximation in real biofilms. This is not a necessary assumption for our model, but it simplifies the equations we consider.

We have reduced the living part of the biofilm, which usually includes many different types of organisms, to one single specie of cyanobacteria, for instance belonging to the Chroococcales order. The four components form a mixture, which can be described as "mixed-state or condition, co-existence of different ingredients or of different groups that mutually diffuse through each other" [START_REF] Rajagopal | Mechanics of mixtures[END_REF]. This approach has been already used by Preziosi et al. to model the formation of vascular tumors [START_REF] Preziosi | On the closure of mass balance models for tumor growth[END_REF].

2.1. Mass balance equations. Since EPS encompasses the cells for this class of cyanobacteria, we can make the hypothesis that live cells, dead cells, and EPS have the same transport velocity, called v S . We denote instead by v L the velocity of liquid phase, and by Γ φ , with (φ = B, D, E, L), the mass exchange rates. Equations expressing the mass balance are

∂ t B + ∇ • (Bv S ) = Γ B , (1a) ∂ t D + ∇ • (Dv S ) = Γ D , (1b) ∂ t E + ∇ • (Ev S ) = Γ E , (1c) ∂ t L + ∇ • (Lv L ) = Γ L . (1d)
We assume the following volume constraint :

(2)

B + D + E + L = 1,
that is, the mixture is saturated and no empty space is left. In addition to the balance of mass of each component, we also have the total conservation of mass of the mixture :

(3)

Γ B + Γ D + Γ E + Γ L = 0,
which states that the mixture is closed, i.e. there is no net production of mass for the whole mixture. Adding the four equations of system (1) and using ( 2) and ( 3) yields

(4) ∇ • ((B + D + E)v S + Lv L ) = ∇ • ((1 -L)v S + Lv L ) = 0,
which means that the divergence of the average hydrodynamic velocity is equal to zero. This can be seen as an average incompressibility.

Observe that when the mass densities are different, but we are still assuming the total mass conservation, we should modify the system by inserting a non null source term in the equation of incompressibility ( 4) with some simple changes in the equation for mass conservation (3).

2.2. Biomass growth rates. We assume that the mass production terms are of the following form :

Γ B = k B LB -k D B, (5a) Γ D = αk D B -k N D, (5b) Γ E = k E LB -εE. (5c)
The term Γ B , i.e. the mass exchange rate for the active bacterial cells, is the difference between a birth term with rate k B L and a death term with rate k D . The birth of new cells at a point depends on the quantity of liquid and of bacteria available in the neighborhood, more precisely bacteria cannot reproduce themselves in absence of liquid, nor of other bacteria. Liquid is therefore a fundamental limiting factor; for simplicity reasons, we choose the rate of the birth term as a linear function of L. The death term in the expression of Γ B gives rise to a creation term in the mass exchange rate for dead cells Γ D , with a coefficient of proportionality α, since a part of the active cells becomes liquid when the cell dies. In Γ D , we also find a natural decay of dead cells with a constant decay rate k N . The production of EPS by active cells is represented by a term with growth rate k E L, while decays occurs at a rate ε. In principle, it should be possible to model the influence on the EPS production of the liquid fraction by using a non-dimensional function f (L) with values between 0 and 1, but here we prefer to restrict our attention to a linear dependence.

We take the mass exchange rate of liquid Γ L in order to enforce condition (3), that is,

Γ L = B ((1 -α)k D -k B L -k E L) + k N D + εE.
The dependence of the growth of biofilms on light will be taken into account through the form of the constant rate k B as explained in the following subsection.

Dependence on light.

Light is a fundamental variable for some types of cyanobacteria (photoautotrophic cyanobacteria), allowing these organisms to photosynthesize inorganic compounds [START_REF] Peeters | The relationship between light intensity and photosynthesis. A simple mathematical model[END_REF]. There also exist a range of temperature and a range of nutrient concentrations necessary to their survival [START_REF] Kuffner | Effects of nitrate, phosphate and iron on the growth of macroalgae and benthic cyanobacteria from Cocos Lagoon, Guam[END_REF]. These environmental factors directly influence the specific growth rates of involved components (B, E, D) and can reduce them when optimal values are not reached (cf. [START_REF] Eilers | A model for the relationship between light-intensity and the rate of photosynthesis in Phytoplankton[END_REF][START_REF] Van Der Grinten | Temperature and lightdependent performance of the cyanobacterium Leptolyngbya foveolarum and the diatom Nitzschia perminuta in mixed biofilms[END_REF][START_REF] Thebault | Comparison between two mathematical formulations of the phytoplankton specific growth rate as a function of light and temperature, in two simulation models (ASTER & YOYO)[END_REF]). We focus here on the dependence on light; the same arguments can be used to take into account the variations of temperature, and the influence of the concentration of nutrients can be computed using classical reaction-diffusion equations. We feel, however, that light is the most important factor to deal with for cyanobacteria and we consider for the moment that temperature and nutrients are not limiting factors; it is foreseen to study their influence in a more complete study [START_REF] Clarelli | Multidimensional model of biofilm growth : stability and influence of environment[END_REF]. Besides, the experimental results about biofilm growth we found in articles like [START_REF] Congestri | Seasonal succession of phototrophic biofilms in an Italian wastewater treatment plant: biovolume, spatial structure and exopolysaccharides[END_REF][START_REF] Pippo | Capsular polysaccharides of cultured phototrophic biofilms[END_REF][START_REF] Stomp | Colorful niches of phototrophic microorganisms shaped by vibrations of the water molecule[END_REF] are mainly obtained with a constant temperature and a sufficient amount of nutrients.

We write the coefficient k B depending on light intensity I as k B (x, y, z,t) = k B0 • g(I(x, y, z,t)), where k B0 is the optimal growth rate, and g(I) ∈ [0, 1] is an efficiency factor. We denote by I 0 (x, y) the light intensity on the upper surface of water, and by I(x, y, z,t) the intensity in the water in the 3D case. We assume that the light intensity is attenuated following the law of photon absorption in the matter [START_REF] Stomp | Colorful niches of phototrophic microorganisms shaped by vibrations of the water molecule[END_REF] as I(x, y, z,t) = I 0 (x, y,t)e -H z µ(s)ds , where µ = µ 0 (1 + µ h (B + E + D)) is the absorption coefficient, µ 0 is the absorption coefficient when the water is clear, and µ h is a coefficient in the biomasses. By experimental observations, an estimate is µ ≈ 0.9 m -1 if the water is turbid, and µ ≈ 0.2 m -1 if the water is clear. We consider finally, following [START_REF] Thebault | Comparison between two mathematical formulations of the phytoplankton specific growth rate as a function of light and temperature, in two simulation models (ASTER & YOYO)[END_REF][START_REF] Eilers | A model for the relationship between light-intensity and the rate of photosynthesis in Phytoplankton[END_REF], an efficiency factor for algae growth of the form [START_REF] Chorin | Numerical simulation of the Navier Stokes equations[END_REF]

g(I) = 2w (1 + β ) Î Î2 + 2β Î + 1
, where Î(x, y, z,t)

= I I opt = I 0 (t) I opt e -H z µ(s,t)ds .
Here I opt is the optimal light intensity, w is the maximum specific growth rate and β is a parameter to be chosen. The function g reaches its maximum value at I opt and vanishes at zero and at infinity. We also assume that the death rate for cyanobacteria k D , the EPS growth rate k E and the decay rate for dead cells k N are independent of I [START_REF] Ni | A Novel Approach to Evaluate the Production Kinetics of Extracellular Polymeric Substances (EPS) by Activated Sludge Using Weighted Nonlinear Least-Squares Analysis[END_REF] in order to keep the model as simple as possible in absence of experimental evidences.

Force balance equations.

In order to find the suitable form of the force balance equations, we follow the works by Preziosi [START_REF] Preziosi | On the closure of mass balance models for tumor growth[END_REF][START_REF] Astanin | Multiphase models of tumor growth[END_REF], where the theory of mixture is used to model tumor growth.

Following [START_REF] Astanin | Multiphase models of tumor growth[END_REF], to determine the change of momentum of one single phase in a reference volume, we have to take into account for four different contributions: the flux momentum through the boundary of the volume, which is the inertial term φ v φ ⊗ v φ ; the internal forces within one phase still through the boundary of the volume, given by the term Tφ , which is the partial stress tensor; the interaction forces among phases in the volume, which give a term of the form mφ ; and finally, the momentum supply related to phase changes in the volume, which gives the term of the form Γ φ v φ . In our case, for simplicity, we neglect the external forces, which can be easily inserted at any time. Therefore, using the divergence theorem, we can write the equation of force balance for the component φ (φ = B, D, E, L) as follows :

(7) ∂ t (φ v φ ) + ∇ • (φ v φ ⊗ v φ ) = ∇ • Tφ + mφ + Γ φ v φ .
From the theory of mixtures [START_REF] Astanin | Multiphase models of tumor growth[END_REF][START_REF] Rajagopal | Mechanics of mixtures[END_REF], we take the standard decomposition of the interaction forces as mφ = m φ + P∇φ , namely, the sum of a friction term m φ between the phases and a term orthogonal to the level curves of the phases P∇φ , where P is the hydrostatic pressure. This term is needed to enforce the saturation condition (2), since otherwise it should be unclear how to verify (4) all the time. In practice, the function P is an unknown of the problem, very much like the pressure in the classical Navier-Stokes equation, and it plays the role of a Lagrange multiplier with respect to the saturation constraint. Following the same idea, we decompose the partial stress tensor as Tφ = -φ PI + φ T φ , where T φ is the excess stress tensor. Equation ( 7) can be rewritten as

(8) ∂ t (φ v φ ) + ∇ • (φ v φ ⊗ v φ ) = m φ -φ ∇P + ∇ • (φ T φ ) + Γ φ v φ .
The total conservation of momentum leads to

(9) ∑ φ ( mφ + Γ φ v φ ) = ∑ φ (m φ + Γ φ v φ ) = 0.
This equation means that the net momentum supply to the mixture due to all the components is equal to zero. If the mixture is closed, we prove that the sum of interaction forces and momentum transfers due to mass exchanges is null. Summing (8) for φ = B, D, E altogether and using ( 2), ( 3) and ( 9), we find

∑ φ =L m φ + Γ φ v φ = -m L -Γ L v L and
using (2) once again, we obtain

∂ t ((1 -L)v S ) + ∇ • ((1 -L)v S ⊗ v S ) = -(1 -L)∇P + ∇ • ( ∑ φ =L φ T φ ) -m L -Γ L v L .
For the liquid phase we have

∂ t (Lv L ) + ∇ • (Lv L ⊗ v L ) = -L∇P + ∇ • (LT L ) + m L + Γ L v L .
To make assumptions on the form of the excess stress tensors, we recall the general form for an elastic fluid, which is

φ T φ = -(Σ(φ ) -φ λ (φ )∇ • v φ )I + φ µ(φ ) ∇v φ + (∇v φ ) T ,
(see for instance [START_REF] Astanin | Multiphase models of tumor growth[END_REF]), formula (39).

In absence of further biological information, we choose for the solid part (φ = B, D, E) to consider a constant stress and we neglect all the shear stress effects. Namely, we choose λ (φ ) = µ(φ ) = 0 and Σ(φ ) = γφ , with γ > 0 to indicate compression. Therefore T φ = -γI and so

∑ φ =L φ T φ = -γ(1 -L)I.
We also assume that T L = 0, namely, that the excess stress tensor is only present in the solid component and that we only find the hydrostatic pressure in the liquid; this means that in absence of bacteria or EPS, the liquid is at rest. This type of assumption is usually adopted in the theory of deformable porous media, where the excess stress tensor T L is neglected in order to get Darcy like laws [START_REF] Rajagopal | Mechanics of mixtures[END_REF].

We finally assume that the interaction forces for the liquid follow the Darcy law [START_REF] Astanin | Multiphase models of tumor growth[END_REF]: this is obtained taking m L proportional to the difference between the relative velocities of fluid and of component, i.e. m L = -M(v Lv S ), where M is an experimental constant.

We can rewrite the equations for the velocities and using (1), (4), we obtain a closed system of equations, including inertial terms :

(10)                          ∂ t B + ∇ • (Bv S ) = B (k B L -k D ) , ∂ t D + ∇ • (Dv S ) = αk D B -k N D, ∂ t E + ∇ • (Ev S ) = k E BL -εE, ∂ t L + ∇ • (Lv L ) = B((1 -α)k D -Lk B -k E L) + k N D + εE, ∂ t ((1 -L)v S ) + ∇ • ((1 -L)v S ⊗ v S ) + (1 -L)∇P = ∇Σ + (M -Γ L )v L -Mv S , ∂ t (Lv L ) + ∇ • (Lv L ⊗ v L ) + L∇P = -(M -Γ L )v L + Mv S , ∇ • ((1 -L)v S + L v L ) = 0.
On the boundary, we impose Neumann conditions for the volume ratios and no-flux boundary conditions for the velocities : [START_REF] Cuzman | Biodiversity of phototrophic biofilms dwelling on monumental fountains[END_REF] ∇B

• n| ∂ Ω = ∇E • n| ∂ Ω = ∇D • n| ∂ Ω = 0, v S • n| ∂ Ω = v L • n| ∂ Ω = 0.
In the following section, we present the numerical scheme we use and the numerical difficulties we had to face to solve numerically this complex system of equations.

NUMERICAL SCHEME

We explain how to solve the complete system (10) by a finite difference method in space and an explicit-implicit method in time in the two-dimensional case. This procedure can be easily adapted to the three-dimensional case.

Consider a square Ω = [0, L] × [0, L] and the discretization grid x α = (α 1 δ x, α 2 δ x), where α 1 , α 2 are integers such that 0 ≤ α 1 , α 2 ≤ N + 1 and δ x is the space step. Denote by δt the time step and by t n = n δt, n ∈ N the discretization times.

Hyperbolic system (10) can be written as

(12) ∂ t W + ∂ x 1 A 1 (W) + ∂ x 2 A 2 (W) = F mass + F f orce + F press ,
where

W =           B D E (1 -L)v S1 (1 -L)v S2 Lv L1 Lv L2           , F mass =           B(k B L -k D ) αk D B -k N D k E BL -εE 0 0 0 0           , ( 13
)
F f orce =           0 0 0 (M -Γ L )v L1 -Mv S1 (M -Γ L )v L2 -Mv S2 -(M -Γ L )v L1 + Mv S1 -(M -Γ L )v L2 + Mv S2           , F press =           0 0 0 -(1 -L)∂ x P -(1 -L)∂ y P -L∂ x P -L∂ y P           , ( 14
)
A 1 (W) =           Bv S1 Dv S1 Ev S1 (1 -L)v S 2 1 + γ(1 -L) (1 -L)v S1 v S2 Lv L 2 1 Lv L1 v L2           , A 2 (W) =           Bv S2 Dv S2 Ev S2 (1 -L)v S1 v S2 (1 -L)v S 2 2 + γ(1 -L) Lv L1 v L2 Lv L 2 2           , (15)
where we used the notations v L = (v L1 , v L2 ) T and v S = (v S1 , v S2 ) T for the twocomponents vectors v S and v L . We split the source term in three different parts in order to make the final scheme we use more readable.

We denote by W n,α the approximation of W at point x α ∈ Ω ⊂ R 2 and time t n . Our final scheme is a modification of the following scheme, which is roughly speaking the Rusanov scheme :

W n+1,α = W n,α - δt 2δ x 2 ∑ j=1 (A j (W n,α j +1 ) -A j (W n,α j -1 )) + λ δt 4δ x 2 ∑ j=1 (W n,α j +1 -2W n,α + W n,α j -1 ) + δt F n,α mass + δt F n,α f orce + δt F n,α press , (16)
where λ is a numerical velocity. This scheme is based on a relaxation scheme for the spatial discretization of the transport part [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. We set the numerical velocity λ as the maximum of the eigenvalues of the Jacobian matrices of the fluxes A 1 and A 2 , namely,

λ = max{2v L1 , v S1 - √ γ, v S1 + √ γ, 2v L2 , v S2 - √ γ, v S2 + √
γ}; space step is fixed, whereas the time step varies and satisfies the stability condition λ δt δ x ≤ 1. We use the Neumann boundary conditions [START_REF] Cuzman | Biodiversity of phototrophic biofilms dwelling on monumental fountains[END_REF] for the components B, D, E, and the no-flux boundary conditions [START_REF] Cuzman | Biodiversity of phototrophic biofilms dwelling on monumental fountains[END_REF] for the velocities v S , v L ; L is computed using (2). This scheme is a standard one; however, considering system (10) and scheme ( 16), we can notice that two difficulties require an original treatment, as explained in details in the two following subsections. The first one deals with the vanishing phases problem, that is, L may be equal to 0 or 1 at some points. Since the scheme computes the quantities (1 -L)v S and Lv L , it is not clear how to compute v S or v L , which are needed in the friction terms. The second one is the computation of the hydrostatic pressure P appearing in the source term F press . We explain in the two following paragraphs how to modify scheme ( 16) in order to face these difficulties.

The vanishing phases problem. Considering scheme [START_REF] Hall-Stoodley | Bacterial Biofilms: from the natural environment to infectious diseases[END_REF] and the definition of W in (13), we note that the computation gives us the two quantities (1 -L n+1 )v S n+1 and L n+1 v L n+1 , whereas we need v S n+1 and v L n+1 in order to compute at time t n+2 the four last components of F n+1,α f orce in [START_REF] Eilers | A model for the relationship between light-intensity and the rate of photosynthesis in Phytoplankton[END_REF] . Therefore, if L n+1 is equal to 1 (resp. 0), we cannot calculate v S n+1 (resp. v L n+1 ). A naive solution would be to take the data in order to avoid vanishing phases, like in most of the works in multiphasic gas dynamic, but those phases are physically relevant in our case, since they correspond respectively to pure liquid or pure biofilm.

Therefore, we propose to treat this problem by using an implicit time discretization for the velocities of the source terms, leaving unchanged the discretization of mass balance equations [START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF]. This leads to replace in [START_REF] Hall-Stoodley | Bacterial Biofilms: from the natural environment to infectious diseases[END_REF] the term F n,α f orce by F n+1,α f orce . This enables us to obtain some values for the velocities even in the particular cases when one of the phases, liquid or solid, vanishes. The only drawback is that we have a system of two non-linear equations to solve; this is always possible under the condition (1 -L n+1, j )L n+1, j + δt(M + (L n+1, j -1)Γ n+1, j L ) = 0, which requires that Γ n+1, j L = M. More precisely, the condition holds true without any restriction if Γ n+1, j L < M, whereas we need to impose δt small enough when Γ n+1, j L > M. The pressure problem. We explain how to compute the term F press . One naive solution would be to find an elliptic equation for P, summing the two force balance equations in [START_REF] Crispim | Cyanobacteria and Biodeterioration of cultural Heritage: A Review[END_REF] and using (4). However, the solution of this equation is not unique and it is not an efficient way to compute the pressure. We will rather use a projection method as the one proposed in [START_REF] Chorin | Numerical simulation of the Navier Stokes equations[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF] to solve Navier-Stokes equations (see also [START_REF] Temam | Theory and numerical analysis[END_REF]). The strategy is to use first a prediction step to compute an estimated velocity and then, in a correction step, to solve an elliptic equation that determines the pressure P, enforcing the divergence constraint (4). The computation of the pressure then enables us to correct the values of the velocities. As a first step, we compute

W n+1/2,α =W n,α - δt 2δ x 2 ∑ j=1 (A j (W n,α j +1 ) -A j (W n,α j -1 )) + λ δt 4δ x 2 ∑ j=1 (W n,α j +1 -2W n,α + W n,α j -1 ) + δt F n,α mass + δt F n+1/2,α f orce .
Note that the source term F f orce is treated implicitly as explained in the previous paragraph. Then, we compute the quantity ∇P n+1 which satisfies the following equation :

δt∆∇P n+1 = ∇∇ • (1 -L n+1 )v S n+1/2 + L n+1 v L n+1/2 ,
with boundary condition ∇P n+1 • n = 0. The final step is to update velocities thanks to the pressure as

v S n+1 = v S n+1/2 -δt∇P n+1 , v L n+1 = v L n+1/2 -δt∇P n+1 .
In the following section, we present numerical simulations in 1D, 2D and 3D. These simulations show the efficiency of our model to reproduce a front propagation for the production of EPS and the growth of the biofilm structure.

NUMERICAL SIMULATIONS

4.1. Parameters estimates. In this section we discuss the values of parameters chosen in our simulations. Some information can be found in [START_REF] Johnson | Adaptive significance of circadian programs in cyanobacteria[END_REF], where a doubling time of the culture of cyanobacteria of once every 10.5 hours is mentioned and [START_REF] Eberl | Mathematical modeling of biofilm[END_REF], where a maximum specific biofilm growth rate is assumed to be 5.88 d -1 , where d = day, while in [START_REF] Thebault | Comparison between two mathematical formulations of the phytoplankton specific growth rate as a function of light and temperature, in two simulation models (ASTER & YOYO)[END_REF] the biofilm growth of the two models presented varies between 1 and 2 d -1 . We choose k B0 = 1 d -1 . According to [START_REF] Eberl | Mathematical modeling of biofilm[END_REF], the decay rate coefficient is about 5% of the growth coefficient and so we chose k D = 0.025k B0 . Finally, the estimate of the EPS growth rate k E depends on the specie of bacteria we study and the environmental conditions (see for example [START_REF] Ni | A Novel Approach to Evaluate the Production Kinetics of Extracellular Polymeric Substances (EPS) by Activated Sludge Using Weighted Nonlinear Least-Squares Analysis[END_REF]). Since we consider the behavior of Chroococcales cyanobacteria, in absence of specific experiments, we assume that the order of magnitude of EPS growth is comparable with the one of cyanobacteria growth.

To estimate the parameters, which are present in the law for the variation of light intensity [START_REF] Chorin | Numerical simulation of the Navier Stokes equations[END_REF], we refer to [START_REF] Eilers | A model for the relationship between light-intensity and the rate of photosynthesis in Phytoplankton[END_REF][START_REF] Thebault | Comparison between two mathematical formulations of the phytoplankton specific growth rate as a function of light and temperature, in two simulation models (ASTER & YOYO)[END_REF]. In [START_REF] Stomp | Colorful niches of phototrophic microorganisms shaped by vibrations of the water molecule[END_REF], a law of absorption of light by water is proposed as well as some estimates for the parameter µ 0 . The parameter µ h is estimated indirectly by some experimental results in [START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF], which link the height of the biofilm and the light transmitted after a path in the biofilm. Experiments in [START_REF] Congestri | Seasonal succession of phototrophic biofilms in an Italian wastewater treatment plant: biovolume, spatial structure and exopolysaccharides[END_REF][START_REF] Pippo | Capsular polysaccharides of cultured phototrophic biofilms[END_REF] suggest the value for the optimal light intensity I opt . The parameter β in ( 6) is estimated thanks to the evolution of the growth as a function of light described in [START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF].

We did not find experimental estimates for the parameter γ in literature. Physically, the value of √ γ plays the role of the sound speed of the mixture, namely, it is the maximum of the modulus of the eigenvalues of the Jacobians of the fluxes when we linearize the system around zero velocities. So, it gives a good estimate of the speed of the front between the solid and the liquid phases. Therefore, according to the experimental results in [START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF], we considered an estimated growth of biofilms of 1 mm in 30 days, which nearly corresponds to a growth of 10 -4 cm in 1 hour. In consequence, for some of the tests we used the value √ γ = 10 -4 3600 [cm/sec] ≈ 2.7 • 10 -8 [cm/sec], which gives γ = 5 • 10 -16 cm 2 /sec 2 . Some numerical estimates of the influence of different values of γ on the growth of the biofilms will be presented later.

Simulations.

In this section we present numerical simulations in one, two and three space dimensions. The one-dimensional system is not a realistic model, but it is useful to test our model on simple cases and to compare it with pre-existing models which have been mostly written only in the 1D context. In this case, the closure law is also more natural

k D 2 • 10 -7 [1/sec] D growth rate [13], k E 12 • 10 -6 [1/sec] EPS growth rate [20] ε 1 • 10 -7 [1/sec] EPS death rate k N 1 • 10 -6 [1/sec] D consumption rate γ 2.5 -10 • 10 -16 cm 2 /sec 2 tensor coefficient α 0.25 [dimensionless] fraction dead cells M 10 -8 [1/sec] tensor coefficient I opt 0.01 µmol • cm -2 sec -1 optimal light intensity [8, 12] I 0 0.01 µmol • cm -2 sec -1 average incident light [8, 12] w 1 [dimensionless] constant light [14, 28] β 0.1 [dimensionless] constant light [34] µ 0 0.002 cm -1 clear water coeff. [25] µ h 6 [dimensionless]
biomasses coeff. [START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF] than in other dimensions; it is also easier to derive an adapted numerical scheme which is a good test case for our modeling and numerical approach.

In the whole section, we use the values of parameters listed in Table 1. The values found for these parameters are only a first approximation; in a further study [START_REF] Clarelli | Multidimensional model of biofilm growth : stability and influence of environment[END_REF], we will analyze the sensitivity of the solutions to all these parameters, showing that the most sensitive one is k B0 which gives a variation of around 10% (below 1.5% for the other parameters) of the volume of biofilms when it varies of 5%.

Finally, in the following subsections we study the influence of light by comparing results obtained with a constant optimal rate for bacteria growth k B = k B0 with a variable k B = k B (I) which depends on light following the law given at section 2.3; in this case, light intensity is taken as a sinusoid of time taking into account the variations between night and day.

4.2.1.

The one-dimensional case. In the one-dimensional case, the single variable space accounts for the height, that is, we consider an homogeneous planar biofilm. We study the evolution of the biofilm after 60 days (see Figure 1). Our domain is the interval L = [0, 1] (cm) and we take spatial step h equal to 0.001cm. We consider as initial data Heaviside functions for cyanobacteria and EPS : B 0 = 0.2 • χ [0,0.007] , E 0 = 0.008 • χ [0,0.007] and the other variables are initially equal to zero. Simulations are done comparing three different values of γ : γ 1 = 10 -16 (cm 2 /sec 2 ), γ 2 = 5 • 10 -16 (cm 2 /sec 2 ) and γ 3 = 10 -15 (cm 2 /sec 2 ) (respectively on the top, in the middle and on the bottom of Figure 1) and two different values of k B : k B equal to the constant k B0 and k B depending on light intensity (respectively, on the left and on the right side of Figure 1).

Let us see how the parameter γ, which is connected to the interaction force among cells and EPS, influences the diffusion of the biomass, that is, the process is accelerated for bigger values of γ : for the biggest value γ = 10 -15 (cm 2 /sec 2 ), the front propagation reaches the highest height (≈ 0.3(cm)), while the heights reached with the values γ = 10 -16 (cm 2 /sec 2 ) and γ = 5 • 10 -16 (cm 2 /sec 2 ) are respectively ≈ 0.1, 0.2. Computing the height of the biofilm is particularly simple in the 1D case, since the sum B + D + E is null for space variable x large enough. A more careful study of the velocity of the front with respect to γ will be performed in the 2D-case.

We can also observe the differences in the cyanobacteria growth using a constant coefficient k B = k B0 or using a coefficient k B = k B (I) which depends on light, which changes according to some periodic function of time. Since k B0 is the optimal growth rate of cyanobacteria, the cyanobacteria growth is higher in the constant case than in the variable light case. We can see in Figure 1 that comparing the curves for a variable light (on the right) and the curves for a constant light (on the left) for the same value of γ, the distribution of liquid is strictly the same. This means that the liquid, which is initially present everywhere, is consumed by the biomass at the same rate in both cases. However, we see a greater growth of EPS in the case of variable light which counterbalances the lower growth of cyanobacteria. 4.2.2. The two-dimensional case. We simulate the biofilm growth in the two dimensional case, during 30 days. In that case, space variables account for length and height, considering that all functions are constant in width. We first study the change of behavior of our model under the variations of light intensity.

Influence of light. We take for B as an initial condition the sum of three Gaussian functions, centered respectively in 0.35, 0.5, and 0.7 (cm) with a total fraction of volume for the cyanobacteria BT 0 = 1.0762 • 10 -5 . Results are presented in Table 2, where function g is defined by [START_REF] Chorin | Numerical simulation of the Navier Stokes equations[END_REF]. We choose an optimal light intensity I opt = 0.01 (µmolcm -2 sec -1 ) and the incident light intensity I 0 defined by [START_REF] Chorin | Numerical simulation of the Navier Stokes equations[END_REF] varies. The effect of light on the volume of biofilms is not so easy to be determined, since the dependence on the light in the model is defined through the form of the coefficients. A plot of the volume fraction of cyanobacteria 2. We observe that the light intensity propagation in the biofilm phase is attenuated by the biofilm mass itself, as expected by the model assumptions. It would be possible to fit the model with this attenuation with appropriate experiments on cyanobacteria growth, measuring the number of photons emitted and the photons which pass through the biofilm with a light sensor under the biofilm. Influence of γ. The parameter γ also has a great and not simply predictable influence on the behavior of the solutions. In Table 3 , we compare two extreme values of γ (maintaining constant all other parameters) for the same initial data as before.

Front velocity. We describe in this paragraph the biofilm front velocity behavior as a function of γ, using as initial condition a single Gaussian function centered in 0.5. In order to define the biofilm front, we consider the ratio between the biofilm volume fraction (that (cm 2 /sec 2 ) (on the top ), γ 2 = 5 • 10 -16 (cm 2 /sec 2 ) (in the middle), γ 3 = 10 -15 (cm 2 /sec 2 ) (on the bottom). On the left, we display the results for an optimal constant rate k B = k B0 and on the right, for a variable rate depending on light periodically in time. is, the sum of B + D + E) and the maximum of this volume fraction at the same time.

We can define two regions, the first one composed of the points where the value of this ratio is less than 0.5% and the second one composed of the points where the value of this ratio is more than 0.5%. The biofilm front will be the boundary between these two regions. We display the movement of the front as a function of time for different values of γ on Figure 3(a). Using a linear fit on the previous curves, we can find the velocity of the front movement, and we plot it as a function of √ γ on the right of Figure 3(b). The curve shows a nice linear behavior. This fact supports our previous claim that the sound speed of the mixture √ γ yields a good estimate of the speed of the front of the solid phase.

Moreover, using a linear fit, we find a relation of the form v f = q 1 √ γ + q 2 , where v f is the front velocity, q 1 = 1.8497 and q 2 = 2.9143 • 10 -8 (cm/sec). This approximation could be compared with experimental data in order to estimate the parameter γ, once the values of other parameters are known.

4.2.3.

The three-dimensional case. In this last subsection, we consider the three dimensional case on a domain

Ω = [0, 1] × [0, 1] × [0, 0.5] (cm 3
) with a sum of Heaviside functions of amplitude 0.02 (cm), as an initial condition for the cyanobacteria volume fraction B and parameter γ = 10 -15 (cm 2 /sec 2 ). Numerical results of the evolution of biofilm with respect to time are displayed in Figure 4 and we observe the formation of an homogeneous layer of the solid components by a quick aggregation, which can be compared with experiments proposed by the Delft's team in [START_REF] Schrauwers | Irregular biofilms are predictable[END_REF]. On the right of Figure 4, we represent simultaneously the evolution of the light intensity I(x, y, z,t), which is attenuated by the biofilm. Although it is not possible to compare our numerical results with experimental data in a quantitative way, we may say that the order of the height of biofilm we find after 30 days (around 1 mm) is in agreement with experiments performed by biologists [START_REF] Pippo | Capsular polysaccharides of cultured phototrophic biofilms[END_REF] on cyanobacteria biofilms under the same conditions of light and temperature than our simulations (around 0.2 mm in [START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF] for slightly different biofilm composition). The qualitative behavior of our biofilm also fits what we know about the evolution of cyanobacteria biofilm. We present finally in Table 4 the computation times needed for the simulations mentioned above using the Matlab software. The space step is equal to 0.001 in the 1D case and to 0.01 or 0.02 in the 2D or 3D cases. In the optimal rate case we mostly used the variable time step, always respecting the CFL condition, which shortens the time needed for the numerical simulations. On the left, we show the evolution of volume fraction of cyanobacteria B and on the right, we present the light intensity variation. On the left, the scale for the volume of cyanobacteria changes since the variation in the volume of biofilms is important, whereas the scale for intensity is fixed. We have proposed a new model of hyperbolic PDEs to describe the formation and growth of a biofilm in a general setting and we have applied it to an environment mostly inspired by the specific case of Chroococcales cyanobacteria. This model is based on the mixture theory with the assumption that bacteria and EPS are moving along together at the same velocity; it naturally reproduces a clear front propagation that has been observed in the development of biofilms, with a finite speed of propagation obtained thanks to inertial terms, which are generally neglected.

A numerical scheme has been produced and tested, where we have treated in a simple way the problems of vanishing phases and of the computation of the pressure term. In order to obtain numerical simulations as realistic as possible, we have estimated the physical parameters from the literature. However a main parameter γ, accounting for the stress function, is missing. We have studied the effect of this parameter, proving that different values give a great variation in the results. A natural and efficient way to calibrate it will be to make a careful comparison with experimental data. We have also studied in details the influence of light which is a crucial environmental factor in the development of phototrophic biofilms. Although the influence of nutrients has been deeply studied in the previous models, it is the first time to our knowledge that light is taken into account for a biofilm model. Numerical simulations show clearly that the variation of light intensity influences the rate of growth of biofilms, which was not obvious from equations since the dependence on the light appears only through the constants.

Finally, the simulations we have performed in the 3D case are clearly comparable with what can be found in the experimental literature [START_REF] Pippo | Capsular polysaccharides of cultured phototrophic biofilms[END_REF][START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF] and give a first validation of our model. For instance, in [START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF], we can find some values of the thickness of the biofilm for different values of the light irradiance. For the moment, no experiments have been performed to compare explicitly the results of this model with experimental data. However, the simplicity of the model would allow us to fit easily experimental results with the parameters. We remark that an important assumption of the experimental data of [START_REF] Zippel | A flow-lane incubator for studying freshwater and marine phototrophic biofilms[END_REF] is the homogeneity of the biofilm, whereas our model shows some more realistic heterogeneities in the biofilm formation.

Thanks to the generality of the model, we should be able in the future to extend it to other applications, such as biofilms with co-existence of various species. Other interesting perspectives would be the study of the influence of signaling between bacteria and also to use the model for preventing infections in a biomedical setting.

FIGURE 1 .

 1 FIGURE 1. Volume fractions of the biofilm components (bacteria in blue, dead bacteria in green, EPS in black and liquid in red) as functions of height (in cm) after 60 days for three different values of γ, γ 1 = 10 -16 (cm 2 /sec 2 ) (on the top ), γ 2 = 5 • 10 -16 (cm 2 /sec 2 ) (in the middle), γ 3 = 10 -15 (cm 2 /sec 2 ) (on the bottom). On the left, we display the

FIGURE 2 .

 2 FIGURE 2. Biofilm evolution, with coefficient γ = 10 -15 , (cm 2 /sec 2 ) and initial condition for B given by the sum of three Gaussian functions centered in 0.35, 0.5, and 0.7 (cm). On the left, we display the volume fraction of cyanobacteria as a function of length and height after 30 days. On the upper right, we present the level curves for function B as a function of length and height and on the lower right, we show the level curves of the light intensity variation (in µmol cm -2 sec -1 ).

FIGURE 3 .

 3 FIGURE 3. On the top, we display the height of the biofilm (in cm) with respect to time (in days) in order to see the biofilm front movement for different values of γ, i.e. γ = (1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5) • 10 -16 (cm 2 /sec 2 ). On the bottom, we plot the front velocity (in cm/s) as a function of √ γ (in cm/s), obtained by a linear fit of the curves on the left.

FIGURE 4 .

 4 FIGURE 4. Numerical simulation of the biofilm at three different times (T = 15, 22, 30 days) with variable rate and with γ = 10 -15 (cm 2 /sec 2 ).On the left, we show the evolution of volume fraction of cyanobacteria B and on the right, we present the light intensity variation. On the left, the scale for the volume of cyanobacteria changes since the variation in the volume of biofilms is important, whereas the scale for intensity is fixed.

TABLE 1

 1

		. A list of (dimensional) parameters
	Param.	Value	Indications	Reference
	k B0	8 • 10 -6 [1/sec]	Cyanob. growth rate [13, 28, 17]

TABLE 2 .

 2 Numerical results in the 2D case : influence of light

	Form of k B Incident light intensity I 0	Total fraction of	Rate of volume growth
			volume of cyanobacteria
	k B0		0.0029	26950 %
	k B0 • g(I) k B0 • g(I) k B0 • g(I)	0.01 0.003 0.001	0.0019 0.0014 2.4481 • 10 -4	17700 % 13060 % 2275 %
	after 30 days is displayed in Figure

TABLE 3 .

 3 Numerical results in the 2D case : influence of γ

	γ	Total fraction of	Rate of volume growth
		volume of cyanobacteria
	1 • 10 -16 1 • 10 -15	0.0027 0.012	25000% 111500%

TABLE 4 .

 4 Computational time for 30 days simulation performed on a standard laptop in various cases

	Dimension Time step Computational time (sec.)
	1D	0.01	194.87
		variable	9.55
	2D	0.02	2386.34
		0.1	684.61
		variable	76.24
	3D	0.02	6616.84
		0.1	1284.44
		variable	690.17

Acknowledgement.

This work has been partially supported by the project "Mathematical problems for the biological damage of monuments" in the CNR-CNRS 2008-2009 agreement, the PRIN project 2008-2009 "Equazioni iperboliche non lineari e fluidodinamica" and by the ANR project MONUMENTALG, ANR-10-JCJC 0103. Most of the computational work was done with the support of CASPUR and Piero Lanucara in the framework of a Standard HPC Grant 2010. Finally, we thank Luigi Preziosi, Antonio Fasano, Piero Tiano, Oana Cuzman, and Francesca Di Pippo for many useful discussions and suggestions. We have a special thank for the late Patrizia Albertano, who was a very enthusiastic guide in this difficult field, and who unfortunately was not able to see the fruits of this collaboration.