C Deygout
email: chloe.deygout@gmail.com.

A Lesne

F Campillo

A Rapaport

A multiscale approach for biofilm modelling: application for growth in a plug flow reactor

Keywords: biofilm, partial differential equations, individual-based model, plug flow reactor, attachment-detachment, advection-diffusion

We propose a new "hybrid" model for the simulation of biofilm growth in a plug flow bioreactor. Our approach consists in combining three scales: a microscopic one for the individual bacteria, a mesoscopic or "coarse-grained" one that homogenises at an intermediate scale the quantities relevant to the attachment/detachment process, and a macroscopic one in terms of substrate concentration. We compare our simulations with the numerical solutions of an existing partial differential equations model at the macroscopic scale.

Introduction

Biofilms are made of complex communities of micro-organisms attached to a surface. They are ubiquitous structures in nature and in industry, where they can have positive roles such as water treatment or negative roles such as material contamination in industrial processes. In the past years, several modelling approaches were used to represent and study biofilms (see Wang and Zhang 2010 for a general review). One approach is a mathematical description, continuous in time and space, where biofilm development is summed up in a few differential equations (see [START_REF] Wanner | Mathematical modeling of biofilms[END_REF][START_REF] Klapper | Mathematical description of microbial biofilms[END_REF] one is interested in following the whole biofilm at a macroscopic scale in terms of its average concentration over a large surface, with the possibility of some mathematical analyses. Analysis can become very difficult when the models try to take into account hydrodynamics (fluid flow) (e.g. [START_REF] Eberl | A new deterministic spatio-temporal continuum model for biofilm development[END_REF][START_REF] Duddu | A two-dimensional continuum model of biofilm growth incorporating fluid flow and shear stress based detachment[END_REF]. The other approach is the use of cellular automata (CA) and individual-based models (IBM). In the former, space is discretised into microscopic-sized boxes and the evolution of the contents of each box is described by simple rules, typically involving the contents of neighbouring cells. In the latter, each individual bacterium is studied and can have its own set of rules to follow (with a combination of deterministic and stochastic rules). Using the IBM approach, one is interested in looking at interactions between the individuals at a microscopic scale (see [START_REF] Hellweger | A bunch of tiny individuals -individualbased modeling for microbes[END_REF]Laspidou et al. 2010 for reviews). Empirical knowledge of the system and its biology can directly be used as the roots of the elementary mechanisms of IBM models. Nevertheless, simulation complexity and runtime can increase very quickly, for example with the number of parameters or that of computer operations. It also means that these models are better suited to study relatively small populations.

In this work, we developed an "hybrid" IBM making use of information from both microscopic and macroscopic scales with coarse-grained descriptions wherever possible (see Figure 1). Indeed, we know that some events at the microscopic scale can appropriately be simplified via mean-field approximations. While there is a continuum between the scales, we chose to apply our approach Fresh growth medium in Effluent out Planktonic cells on a biological system that can be described on three scales spread over this continuum (micro-, meso-and macroscopic). This case study is used to determine the range of validity of our proposed approach. In more complicated situations, where neither a completely microscopic IBM or a PDE system are possible (too many variables or too complicated geometry and boundaries respectively), a mesoscopic approach should prove useful. Also, it should be able to accommodate more easily extensions, such as considering multiple strains in the biofilm.

Biofilm layer

The biological system we describe is a plug-flow reactor, i.e. a long thin tube fed by a laminar flow (Figure 2). It is well suited to study the use of a coarse-grained description for bacterial growth and detachment, in a case where the hydrodynamical aspects are straightforward. The description is pseudo-1D insofar as diffusion that is transverse to the flow is negligible compared to horizontal advection along the flow and the thickness of the biofilm is not explicitly described.

In the coarse-grained IBM, each bacterium has its own coordinate along the continuous space axis, but this space is also divided in discrete boxes, which are the elements used for the coarse-graining approximations. Inside each box along the tube length, we argue that substrate is uniformly distributed, bacteria reaction rates follow a discretised Monod law, regulation of the attachment of planktonic bacteria is well-accounted for by a logistic curve and the attachment of new bacteria is regulated by another function to be described later on. The use of discretised substrate diffusion and advection (with a discretised gradient and Laplace operator), as well as a discretised Monod law for reaction rates are already common [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF]). The IBM spatial boxes have to be big enough for the law of large numbers to hold inside each box but small enough for our homogeneity assumptions regarding substrate concentration to hold as well (Figure 3).

Models

Basic hypotheses

Regarding the physical system, we study a plug flow reactor, i.e. a long and thin tube flow reactor as illustrated on Figure 2 D extends along the z-axis, from z = 0 to z = L. It is fed at z = 0 with growth medium by a laminar flow of fluid in the direction of increasing z and at velocity v which is constant. The external feed contains all nutrients in nonlimiting amounts except one, which is supplied in a constant, growth-limiting concentration s 0 . The flow carries medium, depleted nutrients, bacteria and their by-products out of the reactor at z = L. Substrate s is considered to diffuse with coefficient d s . We assume negligible variation of nutrient concentration transverse to the axial direction of the tube, because D is assumed to be small compared to L.

Regarding the microbial cells, we consider only one bacterial strain. The bacteria can either be suspended in the fluid or attached to the wall. We consider that planktonic bacteria P have an unbiased random motility superimposed on the advection mechanism. It is similar to regular diffusion and characterised by a diffusion coefficient d P . Attached bacteria A are assumed to be immobile and there is a finite carrying capacity of the wall for attachment. Planktonic bacteria might attach themselves to the wall depending on surface availability, while attached bacteria might detach themselves from the wall. Also, when attached bacteria divide, their daughter-cells might either attach themselves or become planktonic depending on surface availability.

The Partial Differential Equations

Our reference PDE model is a system of equations developed by [START_REF] Ballyk | A model of microbial growth in a plug flow reactor with wall attachment[END_REF]. This model is an extension of previous ones describing wall growth in chemostat systems (see Ballyk et al. 2008 for a review of these models) based on the work of [START_REF] Freter | Survival and implantation of escherichia coli in the intestinal tract[END_REF] on the mammalian gut. The purpose of the extension by [START_REF] Ballyk | A model of microbial growth in a plug flow reactor with wall attachment[END_REF] was to account for spatial heterogeneity and material flow. The model accounts for the concentration of substrate s(t, z), the concentration of planktonic bacteria c P (t, z) and the concentration of attached bacteria c A (t, z). Maximum growth rate of j-type bacteria * 1 µ j µ j * 2 1.10 -4 s -1

Half-saturation constant for jtype bacteria * 1

K j K j * 2 1.10 -3 kg.m -3 Maximum amount of attached bacteria concentration c A,∞ [kg.m -2 ] mass m A,∞ [kg] 1.10 -10 kg
Wall fraction occupied by attached bacteria

C A = c A /c A,∞ M A = m A /m A,∞ - Dimensionless coefficient of G(C A ) ǫ - 0.01
Fraction of daughters of attached bacteria finding sites on the wall per surface unit (not per volume unit like the others) as it represents the surface occupied by the attached bacteria. To convert c A (t, z) to the same dimension as c P (t, z), a coefficient called δ is used, which is the ratio of the tube surface to the tube volume (δ = πD πD 2 /4 = 4/D). The equations are the following:

G(C A ) = 1-C A 1+ǫ-C A - - Initial individual mass - m init 1.
∂s ∂t = d s ∂ 2 s ∂z 2 -v ∂s ∂z -c P f P (s)γ -1 -δ c A f A (s)γ -1 ∂c P ∂t = d P ∂ 2 c P ∂z 2 -v ∂c P ∂z + c P (f P (s) -k P ) + δc A f A (s)(1 -G(C A )) -αc P (1 -C A ) + δβc A ∂c A ∂t = c A (f A (s)G(C A ) -k A -β) + αc P (1 -C A )δ -1 (1)
The substrate uptake rates for planktonic and attached bacteria are given by Monod functions:

f j (s) = µ j s K j + s where j = A or P (2)
The fraction of daughter-cells of attached bacteria finding sites on the wall, G(C A ), depends on the amount of surface available for attachment, with (1983) is the following one:

C A = c A /c A,
G(C A ) = 1 -C A 1 + ǫ -C A (3)
where ǫ is typically small. We also have the following Danckwerts boundary conditions [START_REF] Dochain | Dynamical modelling and estimation in wastewater treatment processes[END_REF]:

vs 0 = -d s ∂s ∂z (0, t) + vs(0, t), ∂s ∂z (L, t) = 0 0 = -d P ∂c P ∂z (0, t) + vc P (0, t), ∂c P ∂z (L, t) = 0
and the following initial conditions:

s(z, 0) = s init (z), c P (z, 0) = c P,init (z), c A (z, 0) = c A,init (z), 0 ≤ z ≤ L
These equations lead to two possible steady state regimes: complete washout of the bacteria from the reactor and successful colonisation of the reactor by the bacteria. Details about these steady states can be found in [START_REF] Ballyk | A model of microbial growth in a plug flow reactor with wall attachment[END_REF]. In order to study the system more closely, a numerical approach was used to solve it. We used a similar approach as that used by [START_REF] Ballyk | A model of microbial growth in a plug flow reactor with wall attachment[END_REF] with a centered second-order finite difference scheme for diffusion, a backward first-order finite difference scheme for advection and the semi-implicit Crank-Nicholson method for time.

The microscopic rules

We are interested in describing attachment processes of our system and understanding the impact and relevance of our description choices. Indeed, studying the results of the PDE system shows that the choice of the G function for the attachment of daughter-cells of attached bacteria and the attachment/detachment rates have an important impact in defining the location, shape and length of the biofilm. There are two individual-level mechanisms that involve bacterial attachment: attachment of planktonic bacteria and attachment of daughter-cells of attached bacteria. The space available for attachment is assumed to be limited. The most basic description of these processes on a microscopic scale would be the following: any arriving bacterium attaches itself with probability α if this bacterium is above or next to an available position, else the bacterium becomes or remains planktonic. This mechanism is the most "microscopic" one we consider, where the complex physical and biological attachment process itself is summed up inside the α parameter. In order to simplify neighbourhood references, it is possible to homogenise such a mechanism at the mesoscopic scale. Thus, there is no more need to follow exactly the trajectory of each bacterium and check whether the position below it is available or not. Let us consider there are X max individual positions on the surface in a given neighbourhood, the elementary volume, and already X individuals occupying some of these positions. A first possibility would be: any arriving bacterium attaches itself with probability α if there are available positions in the considered space (X < X max ), else the bacterium becomes or remains planktonic. This "homogenised threshold" mechanism (HT) is the most simple description at the mesoscopic level of an elementary volume, without density-dependence effects. This also means that, using this mechanism assumes bacteria can attach relatively easily, as long as there is some positions available for attachment. Another more complex description would be: any arriving bacterium attaches itself with probability α(1 -X/X max ), which means that the bacterium becomes or remains planktonic with probability αX/X max . This "homogenised probability" mechanism (HP) is equivalent to the microscopic mechanism, with regard to the bacteria macroscopic behaviour.

It is interesting to notice that in their equations, [START_REF] Ballyk | A model of microbial growth in a plug flow reactor with wall attachment[END_REF] use different mechanisms for the attachment of planktonic bacteria and for the attachment of daughter-cells of attached bacteria. The former is a macroscopic version of the HP mechanism while the latter is of the HT mechanism. These choices are not discussed much in this work nor in their references [START_REF] Ballyk | A model of microbial growth in a plug flow reactor with wall attachment[END_REF]. Thus, we have chosen to explore these different mesoscopic mechanisms in our coarse-grained IBM and their impact on model results. Unless specified otherwise, the IBM runs by default with the HT mechanism.

The coarse-grained Individual-Based Model

The model description follows the ODD (Overview, Design concepts, Details)

protocol for describing individual-and agent-based models [START_REF] Grimm | A standard protocol for describing individual-based and agent-based models[END_REF][START_REF] Grimm | Railsback. The ODD protocol: A review and first update[END_REF]. "Overview" is covered in Sections 2.4.1 to 2.4.3, "Design concepts" in Section 2.4.4 and "Details" in Sections 2.4.5 to 2.4.7.

Purpose

The purpose of this IBM is to represent the biological system described in Section 2.1 while taking into account information from both the macroscopic PDE and a microscopic IBM. We use a "simple" case study to determine the range of validity of this approach.

Entities, macroscopic state variables and scales

The entities in this model are bacteria. Their microscopic attributes are their position, mass and attachment status. Our main state variables are the number of bacteria in the tube, both attached or detached, as well as the concentration of the substrate. The substrate is not described at a molecular level but its concentration is an attribute of each grid box. The time scale of our simulations is about a few days in duration with a few seconds for the time step. Simulation duration is limited to the time it takes to reach a quasi-steady-state in the system. The time step is constrained by the discretisation scheme used to model the substrate diffusion and advection. The space scale is a few centimetres for the tube length with about 50 micrometers for the space step.

Process overview and scheduling

The first process is substrate reaction-diffusion. The substrate in each grid box diffuses horizontally into the two neighbour boxes (the model being pseudo-1D, we assume homogeneity in the directions transverse to the tube axis) and is advected along the flow into the following box (discretised Laplacian). Bacteria in the box eat one after the other in a random order. The speed with which bacteria eat depends on the local substrate concentration and is given by the Monod law (Eq. 2) as long as substrate is available in the box. The amount eaten by the bacteria increases their mass by that same amount multiplied by the yield factor γ. Substrate concentration in the boxes is updated only at the end of these steps. Despite the artificiality of this sequence, discrepancies should be small because the time step remains very small.

The second process is bacteria division. Bacteria that have a mass higher that the division threshold m div divide into two bacteria of total mass that of the mother bacteria (see 2.4.7 for details).

The third process is bacteria attachment/detachment. Finally, the detached bacteria are advected by the flow and have a random horizontal motility.

Design concepts

Basic principles. This work is based on that of [START_REF] Mabrouk | Bacteria can form interconnected microcolonies when a self-excreted product reduces their surface motility: evidence from individual-based model simulations[END_REF]; [START_REF] Mabrouk | Analyzing individual-based models of microbial systems[END_REF], itself based on similar assumptions as that of [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF]; [START_REF] Kreft | Individual-based modelling of biofilms[END_REF]. Main differences regard the presence of a flow and taking the border into account.

Adaptation, objectives, learning, prediction. Given the small time scale of our simulations, the bacteria in this model do not have adaptive traits and hence do not attempt to meet objectives. They have no learning abilities and do not attempt to predict their future environment.

Sensing. The internal variable taken into account by the bacteria is their own mass (which determines division time). The environmental variables taken into account are local substrate concentration (only that of the box in which the bacterium is present), local bacteria abundance (only that of the box in which the bacterium is present) and flow velocity (constant throughout the tube).

Interaction. The interactions between the bacteria are only indirect ones, through competition for resources: substrate and attachment sites.

Stochasticity. Several events in the model are described with probability rules. When a bacterium divides, the mass and location of the daughter bacteria have a random component (see above). Attachment and detachment are also random (see above).

Collectives. When a box reaches its carrying capacity (total attached mass in the box superior or equal to m A,∞ ), the bacteria form a biofilm on the wall. This can be considered as a collective. Once this collective is formed, new bacteria cannot attach themselves in that box. Nevertheless, this collective is not permanent, and if enough detachment occurs, the mass may again become lower than the threshold m A,∞ .

Observation. Data are collected as the last process in the schedule, for each box. We gather information on substrate concentration, bacteria mass and numbers (both attached or detached), attachment and detachment events.

Initialisation

Before the simulation starts, a number of bacteria are placed in the boxes (location picked from a random uniform distribution), a certain number of which are attached (N A,init ), the others detached (N P,init ). Each bacterium has an initial mass of m init . Substrate is present throughout the tube at a given concentration s init .

Input data

The input data, such as bacteria initial position, are chosen randomly. Nevertheless, if one wanted to specifically reproduce a given part of a larger tube, initial conditions of the model could be forced to resemble any point of the PDE system representing the whole tube.

Submodels

Substrate diffusion, advection and reaction. All these processes are mingled in the model and space is considered sequentially along the flow direction.

1. First, we consider the substrate reaching equilibrium where the flow arrives. Thus, the concentration of substrate in box i, s t,i , increases with substrate arriving by diffusion and advection from the upflow box i -1 and decreases from substrate leaving by diffusion into box i -1.

2. Then, we determine the amount of substrate which is eaten by the bacteria inside box i (reaction) based on the amount determined in step 1.

3. Finally, we consider the substrate reaching equilibrium where the flows leaves the box. The concentration of substrate in box i at t + dt, s t+dt,i , is the amount determined in step 2 added with substrate arriving by diffusion from the downflow box i + 1 and once removed substrate leaving by diffusion and advection into box i + 1.

Bacterial biology and movement.

Growth Each bacterium consumes a given amount of substrate which is determined by the Monod growth kinetics and its own mass m. We assume that bacterial growth is directly equal to this amount consumed, multiplied by the conversion yield factor γ. At time t and in box i, the increase in mass dm is the following:

dm = γm µ j s t,i K j + s t,i
where j is either A or P for attached or planktonic bacteria respectively (see Table 1 for more details).

Division Each bacterium that has a mass higher than the division threshold m div divides into two bacteria of total mass that of the mother bacterium.

Each of the two daughter-bacteria receives half of the mother mass plus or minus a random amount (taken from a uniform distribution between 0 and 25% of the mother mass). One of the daughter bacteria retains the mother location and attachment status, and the other bacterium has the same location plus or minus a random amount (taken from a uniform distribution between 0 and mother diameter). That second daughter bacterium has the same attachment status as its mother by default. If that status is attached, depending on the mechanism used for attachment (see Section 2.3), then the new bacterium might become detached. A daughter bacterium which location is outside the tube is be considered washed out.

Attachment and detachment Each of the detached bacteria has a probability involving α times dt and m A,∞ to attach itself on the wall. The exact mechanism can be varied according to the description in Section 2.3. Each of the attached bacteria has a probability β times dt to detach itself.

Movement The detached bacteria are advected by the flow and have a random horizontal motility modelled as a Brownian motion process with an apparent diffusion factor d P . The distance crossed over one time step thus follows a Gaussian law of variance twice the diffusion-like factor d P multiplied by the time step (Mabrouk 2010).

Numerical experiments

We have run a large number of numerical experiments for both the PDE system and the IBM. Parameter values are those indicated in Table 1 unless mentioned otherwise. In the case of the PDE system, one run gives the deterministic output for a given set of parameters. But in the case of the IBM, several runs were done with different random seeds in order to obtain means of the output variable to represent the average behaviour of the model. Because the stochastic aspects of the IBM are small-scale events, a relatively low number of repetitions was sufficient to average out their effects (a dozen runs). The results we study are: the profiles along the tube for all three variables (in terms of mass, for comparison between the two models), the biofilm length defined as the location along the tube where the difference between two following positions is the largest over the whole profile, and the sum over the whole tube of the attached biomass.

Results

General description of steady state profiles

In both models, we find similar profiles at the steady state for the default values (Figure 4). These profiles indicate an important development of the biofilm in the leftmost side of the tube where the nutrient flow enters. Parallel to this biofilm development which reaches the carrying capacity, there is an important decrease in substrate concentration. Once a threshold substrate concentration is crossed, there is not enough substrate to allow the growth of a thick biofilm and the amount of attached bacteria reaches a value close to one bacterium along the remainder of the tube. Meanwhile, the amount of detached bacteria increases along the whole tube, these bacteria being released by the thick biofilm of the leftmost part of the tube.

There are nevertheless differences between the two models. First of all, in the right-most part of the tube, there are few attached bacteria and the two models represent this situation differently. The IBM gives a presence-absence tube length masses (in kg, log-scale)

Figure 4: General profile along the tube observed at the steady state for the PDE numerical resolution (dashed lines) and the average over 12 IBM runs (solid lines) with default values as shown in Table 1. The mass of attached bacteria m A is represented in blue, the mass of planktonic bacteria m P in green and the mass of substrate m S in red. The horizontal axis is the length of tube.

information: there is no point in case no attached bacteria are present due to the log-scale of Figure 4. Meanwhile, the PDE gives a continuous expression which in the right-most part of the graph corresponds to about a tenth of an individual bacterium. Nevertheless, both cases give similar values. Second, the attached bacteria in the thick biofilm are slightly more developed in the IBM than in the PDE, which implies a slightly shorter biofilm length. The difference in biofilm length is exaggerated by the log-scale but it remains quite small (about 1 "box" of difference) and the difference in overall attached biomass is about 2%. Finally, there is also the effect of fluctuations in the IBM to be considered. This variability does not reflect individual variance because there has already been some homogenisation but it could be thought of as a "quality" assessment of the coarse-grained IBM. Overall, the variance along the profile is relatively low, in particular for the substrate concentration and the planktonic bacteria. Indeed, these two elements are affected by diffusion processes which smooth their distribution along the tube and there is less variability between runs. For the attached biomass, it is quite low over the thick biofilm length but increases afterwards, in particular over the switch from thick biofilm to almost no attached bacteria (data not shown).

Parallel between discretisation to solve PDE and discretised events in IBM

In both cases, whether the model was conceived in a continuous or a discrete manner, the resulting simulations relied on a discretised code. Thus, in both cases, we can compare the impact of the number of spatial boxes making up the whole tube. By keeping the tube length constant and increasing the number of boxes, the size and volume of each box is reduced. Using a large number of boxes increases the similarity of the results with those that might be obtained with a continuous representation. But it also increases simulation run time and it reduces the number of individuals or biomass present in each box which may lead to observing quantities smaller than one individual in the case of the PDE system. In this work, we find that the numerical experiments are relatively resilient to the number of boxes and thus to box size as well (Figure 5). By using a small number of boxes, there is a slight tendency to overestimate the biofilm length in both models, but the difference remains well under 1%.

Sensitivity analysis

Overall, the same effects are observed on PDE and IBM simulations. Differences in the main results (length of biofilm, mass of attached bacteria) between the two models remained well under 5%. The parameters tested here are those which have a major impact on the results. An increase in s 0 leads to an increase in the length of biofilm formed (Figure 6). An increase in v leads to a similar effect. Indeed, increasing v leads to an increased amount of substrate available for the bacteria.

An increase in both µ j or both K j leads first to an increase in the length of biofilm formed and after reaching a peak value it leads to a decrease (Figure 7). Indeed, for low values of µ j , bacteria are washed out while for high values of µ j , the bacteria are very efficient in eating the substrate and only a short length of biofilm can develop before the tube is emptied of substrate. Thus, it is for intermediate values of µ j that a maximum length of biofilm is observed. Similar results are observed for K j but with high values of K j leading to a wash-out and low values of K j leading to a short and efficient biofilm.

An increase in m A,∞ leads to a decrease in the length of biofilm which is formed (Figure 8). Indeed, since the biofilm is thicker, the substrate is completely eaten earlier along the tube and prevents further biofilm development.

An increase in α has very little effect on the development of the main biofilm in the leftmost part of the tube. Nevertheless, it leads to an increase in the development of the biofilm in the rest of the tube (Figure 9). Bacteria in this part of the biofilm do not grow because there is very little substrate left but the biofilm is built via attachment of planktonic bacteria flowing by. An increase in β leads to a decrease in the length of biofilm formed (Figure 9). Indeed, for large values of β, the amount of substrate needed for biofilm development (including compensation for biofilm detachment) is higher than for smaller values of β. For a large enough value of β, there is wash-out.

Impact of microscopic attachment mechanisms

Different attachment mechanisms were considered (Figure 10). The PDE mechanism of the G function and the "homogenised threshold" (HT) mechanism lead to similar macroscopic behaviours, with an abrupt limit to attachment possibilities for bacteria and thus an abrupt distinction between a part of the tube with a thick biofilm and the other part with almost no attached bacteria.

Whichever mechanism is followed for the attachment of planktonic bacteria, the results are similar because this attachment is mostly important in the part of the tube with little attached bacteria and thus for values where all mechanisms behave similarly. On the other hand, whether the attachment of daughter bacteria follows the HT or the HP mechanism has a much more important impact In the IBM, the PDE mechanism of the G function and the "homogenised threshold" (HT) mechanism give similar results whether used for the attachment of planktonic or daughter bacteria and their results are thus both represented by the dotted lines. The dashed lines represent the use of the "homogenised threshold" (HT) mechanism for the attachment of daughter bacteria and of the "homogenised probability" (HP) for the attachment of planktonic bacteria in the IBM. The dash-dot lines represent the use of the "homogenised probability" (HP) mechanism for both kinds of attachment in the IBM. The mass of attached bacteria m A is represented in blue, the mass of planktonic bacteria m P in green and the mass of substrate m S in red. The horizontal axis is the length of tube.

on the overall profile along the tube. Indeed, the "homogenised probability" (HP) mechanism leads to a much smoother transition between the part of the tube with a thick biofilm and the other part with almost no attached bacteria. It remains to be seen which may be the more realistic approach.

Discussion

Modelling choices

First of all, our model represents a 3-dimensional system in a pseudo-1dimensional manner. We take only the coordinate along the tube axis (z-axis) into account, assuming that at a given z-coordinate substrate concentration is homogeneous inside the tube and that the precise location of bacteria is much less important than whether they are attached or planktonic. In general, it has been observed that minimal models offer an interesting trade-off between realism and generality.

In this work we have built an IBM with the aim to represent a given biological system and then justify choices made in the pre-existing PDE system. In order to build the IBM, a number of elements could not be derived straightforwardly. For example, all kinds of stochastic events had to be added, in particular regarding the division mechanism. We chose to use individual mass as a nonstochastic indicator of division time, but the distribution of this mass between the two daughter bacteria had a stochastic component. This simple mechanism proved sufficient to prevent any synchronisation among the bacteria and thus more complex ones were not used in the IBM. An important part of this work relied on the study of attachment patterns and various mechanisms were proposed both for attachment of planktonic bacteria and daughters of attached bacteria. Following the equation system, we studied both kinds of attachment separately. Nevertheless, one might wonder whether it is relevant to distinguish between them. Thus, we might suggest using the same mechanism for both, keeping in mind that in this particular system the precise mechanism behind the attachment of planktonic bacteria has a much minor effect than that of the daughters of attached bacteria.

Perspectives

Using an IBM approach to describe this plug-flow system has several advantages. First of all, it allows us to use "biological" rules, i.e. rules that can intuitively be suggested (and tested) by biologists based on their in-depth knowledge of the real system. This makes the system more amenable to interactions with experiments and also to changes in the choice of mechanisms.

This IBM may thus lead to improvement over the original model proposed by [START_REF] Freter | Survival and implantation of escherichia coli in the intestinal tract[END_REF] and its variations [START_REF] Ballyk | A model of microbial growth in a plug flow reactor with wall attachment[END_REF] by allowing us to test mechanisms that could not be easily formulated or even solved analytically. One possibility would be to explore the biofilm dynamics when considering different bacterial strains competing inside the tube for space and/or substrate resources. Also, we could consider taking into account the (transversally) inhomogeneous velocity profile through a shear rate, generating a shear force on attached bacteria which would promote their detachment (with an effect on the β parameter) and limit the thickness of biofilm (effect on the c A,∞ carrying capacity parameter).

Regarding the density issue (thick biofilm in one part of the tube, few attached bacteria in the other part), we could consider an association of the different models. Indeed, currently we have compared the PDE and IBM models but a prospect could be to use them sequentially with the PDE at the start of the tube where the population is larger and a transition to the IBM for the rest of the tube where the population is smaller.

Finally, with this coarse-grained IBM, we open discussions on top-down interactions (macroscopic constraints on individual bacteria or mesoscopic boxes, feedbacks) and on bottom-up interactions (emergent collective features).

Conclusion

We have proposed a new multiscale model of biofilm growth in a plug flow reactor, with reasonable runtimes for the simulations while remaining close to microscopic processes and knowledge. In particular, our approach focuses on the hypotheses regarding the attachment/detachment process that is at the core of biofilm growth. The comparison with the PDE model has led to similar steadystate concentration profiles at the macroscopic level, justifying qualitatively the attachment/detachment terms introduced at the mesoscopic scale, although some differences could be observed. An IBM approach such as the one we use provides additional information on second-order terms or variance around the average profile that have been compared with the PDE model. Our framework is suitable for extensions to multi-species biofilm or more complex interactions between individuals.

Figure 1 :

 1 Figure 1: Double modelling diagram.

Figure 2 :

 2 Figure 2: Flow reactor with biofilm.

Figure 3 :

 3 Figure 3: Definition of the different scales.

  ∞ being the fraction of surface already occupied. The main biological requirement on this function is to decrease when C A increases and reach zero when C A is null. The particular function used since the work of Freter et al.

Figure 5 :

 5 Figure 5: Effect of the number of boxes on the steady state profile for (a) the PDE numerical resolution and (b) the average over ... IBM runs. The mass of attached bacteria m A is represented in blue, the mass of planktonic bacteria m P in green and the mass of substrate m S in red. The horizontal axis is the box coordinate along the tube divided by the number of boxes in order to compare between cases.

Figure 6 :

 6 Figure6: Effect of the input concentration of substrate s 0 on the steady state profile for (a) the PDE numerical resolution and (b) the average over 12 IBM runs. The mass of attached bacteria m A is represented in blue, the mass of planktonic bacteria m P in green and the mass of substrate m S in red. The horizontal axis is the length of tube.

Figure 7 :Figure 8 :

 78 Figure7: Effect of (i) both µ j or (ii) both K j on the steady state profile for (a) the PDE numerical resolution and (b) the average over 12 IBM runs. The mass of attached bacteria m A is represented in blue, the mass of planktonic bacteria m P in green and the mass of substrate m S in red. The horizontal axis is the length of tube.

Figure 9 :

 9 Figure9: Effect of (i) α and (ii) β on the steady state profile for (a) the PDE numerical resolution and (b) the average over 12 IBM runs. The mass of attached bacteria m A is represented in blue, the mass of planktonic bacteria m P in green and the mass of substrate m S in red. The horizontal axis is the length of tube.

  The latter concentration is measured as a weight

	Description		PDE symbol	IBM symbol	Default
					value
	Space		continuous posi-	box number i,	200 boxes
			tion z [m]	space step dz [m]
	Time [s]		continuous time	discrete time t,	-
			t	time step dt
	Substrate	concentration	s(t, z)	s t,i	-
	[kg.m -3 ]
	Amount of planktonic bacteria	concentration	number N P,t,i or	-
			c P (t, z) [kg.m -3 ]	mass m P,t,i [kg]
	Amount of attached bacteria	concentration	number N A,t,i or	-
			c A (t, z) [kg.m -2 ]	mass m A,t,i [kg]
	Initial substrate concentration	s init	s init	1.10 -3 kg.m -3
	Initial amount of j-type bacteria	concentration	number N j,init	1000 ind.
	* 1		c j,init [kg.m -2 ]
	Length of tube		L	L	0.01 m
	Diameter of tube		D	D	1.10 -3 m
	Flow velocity		v	v	1.10 -6 m.s -1
	Conversion yield		γ	γ	1
	Ratio of tube circumference to	δ	-	4.10 3 m -1
	cross-sectional area
	Input substrate concentration	s 0	s 0	1.10 -3 kg.m -3
	Attachment rate		α	α * 2	1.10 -8 s -1
	Detachment rate		β	β * 2	1.10 -5 s -1
	Diffusion rate of substrate	d s	d s * 2	1.10 -14 m 2 .s -1
	Diffusion rate of planktonic bac-	d P	d P * 2	1.10 -13 m 2 .s -1
	teria

Table 1 :

 1 Description of the symbols used in both models, with units shown in brackets. Footnotes: * 1 j = P for planktonic bacteria or A for attached bacteria; here P and A parameters have the same values. * 2 Same coefficient in both models but may be used slightly differently to account for discretisation.

	10 -15 kg

Acknowledgments

This work has been achieved within the DISCO ("multi-scale modelling bio-DIversity Structure COupling in biofilms") project of the French Agency of Research SYSCOMM program (ANR AAP215-SYSCOMM-2009), that has financed CD's postdoctoral fellowship.


Marcelo Freitas Caetano

Xavier Rodet

Improved Estimation of the Amplitude Envelope of Time Domain Signals Using

Keywords: Amplitude envelope estimation, true envelope, cepstral smoothing, musical instrument sound, speech

whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The estimation of the amplitude envelope of a time-domain signal (or waveform) is a classical engineering problem that arises in amplitude demodulation [START_REF] Rodet | A Comparison of the Energy Operator and the Hilbert Transform Approach to Signal and Speech Demodulation[END_REF], onset detection [2], and attack time estimation [START_REF] Luce | Durations of Attack Transients of Nonpercussive Orchestral Instruments[END_REF], temporal modeling of sounds such as automatic transcription [START_REF] Schloss | On the Automatic Transcription of Percussive Music -From Acoustic Signal to High-Level Analysis[END_REF] and temporal evolution of musical instrument sounds [START_REF] Burred | Automatic Segmentation of the Temporal Evolution of Isolated Acoustic Musical Instrument Sounds Using Spectro-Temporal Cues[END_REF]. The amplitude envelope is a perceptually important feature of musical instrument sounds and speech. It has been shown to be correlated to the percussiveness of musical instruments sounds [START_REF] Mckinney | Features for Audio Classification: Percussiveness of Sounds[END_REF], to speech intelligibility [START_REF] Drullman | Temporal Envelope and Fine Structure Cues for Speech Intelligibility[END_REF] and even to affect pitch perception [START_REF] Hartmann | The Effect of Amplitude Envelope on the Pitch of Sine Wave Tones[END_REF]. The classical approach to the estimation of the amplitude envelope of a time-domain signal is the technique known as envelope follower [START_REF] Rodet | A Comparison of the Energy Operator and the Hilbert Transform Approach to Signal and Speech Demodulation[END_REF], [START_REF] Schloss | On the Automatic Transcription of Percussive Music -From Acoustic Signal to High-Level Analysis[END_REF], [START_REF] Mckinney | Features for Audio Classification: Percussiveness of Sounds[END_REF], which consists basically of rectifying the waveform and then low-pass filtering it. This can be easily implemented in either the analog or digital domains. Another very popular solution is to calculate the instantaneous root mean square (RMS) value of the waveform through a sliding window with finite support [START_REF] Hajda | A New Model for Segmenting the Envelope of Musical Signals: The relative Salience of Steady State versus Attack, Revisited[END_REF]. More recently, some authors have proposed other techniques to obtain a more reliable estimation of the amplitude envelope of waveforms. An early attempt [START_REF] Schloss | On the Automatic Transcription of Percussive Music -From Acoustic Signal to High-Level Analysis[END_REF] consisted of a piece-wise linear approximation of the waveform. The amplitude envelope is created by finding and connecting the peaks of the waveform in a window that moves through the data. Jensen [START_REF] Jensen | Envelope Model of Isolated Musical Sounds[END_REF] proposed a method that fits curve shape approximations to model the amplitude envelope of the partials of an additive model of instrument sounds and later Skowronek [START_REF] Mckinney | Features for Audio Classification: Percussiveness of Sounds[END_REF] applied it to approximate the global amplitude envelope. However, advances in spectral envelope estimation techniques such as linear prediction [START_REF] Makhoul | Linear prediction: A Tutorial Review[END_REF] and cepstral smoothing [START_REF] Villavicencio | On Cepstral and All-Pole Based Spectral Envelope Modeling with Unknown Model Order[END_REF] are promising candidates to obtain an accurate estimation of the amplitude envelope when applied in the time-domain [START_REF] Burred | Automatic Segmentation of the Temporal Evolution of Isolated Acoustic Musical Instrument Sounds Using Spectro-Temporal Cues[END_REF]. In this work we propose to apply the true envelope estimation technique in the time domain to accurately estimate the amplitude envelope of waveforms such as isolated musical instrument sounds or speech. We will show that the true amplitude envelope renders estimations that respond well to sudden changes in amplitude while remaining smooth during more stable regions of the waveform.

In the next section we present the classical amplitude envelope estimation techniques, followed by the recently proposed use of linear prediction. Next, we introduce the true envelope cepstral smoothing algorithm and we explain how we apply it in the time domain. Finally, we present and evaluate the results, comparing our proposed approach with the others.

CLASSICAL AMPLITUDE ENVELOPE ESTIMATION

The classical amplitude envelope estimation techniques explained in this section are low-pass filtering (LPF), root-mean square (RMS), and analytic signal. We will also review a recent proposal to the use of linear prediction [START_REF] Ellis | Frequency-Domain Linear Prediction for Temporal Features[END_REF], dubbed frequency-domain linear prediction (FDLP).

Low-Pass Filtering (LPF)

Low-pass filtering is the most straightforward way of obtaining a smooth signal that follows the amplitude evolution of the original waveform. It is based on a classical amplitude demodulation envelope follower technique [2], that low-pass filters a half-wave (hwr) or full-wave rectified (fwr) version of an amplitude modulated (AM) signal. The principle of amplitude modulation (AM) is that the amplitude changes of the signal carry the information we seek. There are many possible filter designs with different characteristics and the choice affects the quality of the final envelope. For instance, Jensen [START_REF] Jensen | Envelope Model of Isolated Musical Sounds[END_REF] proposes to convolve the waveform with a Gaussian window function, resulting in a suboptimal estimation. Also, the cut-off frequency of the filter has a major impact on the result. High cut-off frequencies will likely produce an amplitude envelope with ripples and very low cut-off frequencies are less responsive to sudden amplitude changes.

Root-Mean Square (RMS) Energy

The RMS value is perhaps the most popular [START_REF] Hajda | A New Model for Segmenting the Envelope of Musical Signals: The relative Salience of Steady State versus Attack, Revisited[END_REF] method for estimating the temporal evolution of the signal energy because it can be easily adapted to obtain an estimate of the amplitude envelope by simply applying it with a sliding window, as shown in equation (1) ( ) ( ) ( )

t x t w T t RMS T i i i ∑ = = 1 2 1 ( 1
)
where x i (t) is the i th sample of the signal centered around t as seen through the window w i (t), t is the number of samples the analysis window moves, and T is the window length. Usually a rectangular window is used, but other choices are also possible [START_REF] Jensen | Envelope Model of Isolated Musical Sounds[END_REF]. The RMS calculation is a special case of the generalized mean with exponent p=2 and as such, also functions as a sort of moving average, lowpass filter that smoothes out the signal. The analysis step t imposes a trade-off between the temporal sampling rate of the envelope and how much information it represents. Small values of t react sooner to sudden changes in amplitude, while presenting ripple in more steady regions and larger values smooth out the ripples but tend to lag behind abrupt energy changes.

Analytic Signal

The Hilbert transform is part of a signal processing technique for amplitude demodulation [START_REF] Rodet | A Comparison of the Energy Operator and the Hilbert Transform Approach to Signal and Speech Demodulation[END_REF]. The Hilbert transform of a signal x(t) is defined as

( ) ( ) t t x t x π 1 ˆ * = (2)
where * stands for convolution. Using equation ( 2), we can define the analytic signal z(t) as

( ) ( ) ( ) ( ) ( ) [ ] t j t r t x j t x t z θ exp ˆ= + = (3).
The analytic signal is useful for envelope detection since its modulus r(t) and time derivative of the phase θ(t) can serve as estimates for the amplitude envelope and instantaneous frequency of x(t) under certain conditions. Notably, if the Hilbert transform of x(t) is equal to its quadrature signal [START_REF] Rodet | A Comparison of the Energy Operator and the Hilbert Transform Approach to Signal and Speech Demodulation[END_REF], then the estimates are equal to the actual information signals [START_REF] Rodet | A Comparison of the Energy Operator and the Hilbert Transform Approach to Signal and Speech Demodulation[END_REF]. Synthetic (i.e., AM) signals can be constructed to have this property, but there is no reason to expect that acoustic musical instrument sounds or speech also present it. A more realistic condition is verified when we are dealing with narrowband signals [START_REF] Hartmann | Signals, Sound, and Sensation[END_REF], which is rarely the case for musical instrument sounds and speech. The analytic signal can be effectively used to extract the amplitude envelope of individual partials if applied to each frequency bin of the STFT, but when applied to the whole signal it is equivalent to trying to demodulate several AM signals at the same time, so we use it as half-wave rectifier in this work.

Frequency-Domain Linear Prediction (FDLP)

Traditional linear prediction [START_REF] Makhoul | Linear prediction: A Tutorial Review[END_REF] estimates the spectral envelope from the time-domain signal. The idea behind FDLP [START_REF] Ellis | Frequency-Domain Linear Prediction for Temporal Features[END_REF] is to exploit time-frequency duality to extract the temporal amplitude envelope by applying linear prediction to a spectral representation. In particular, the used spectral representation is the discrete cosine transform (DCT) given by equation ( 4), since it is real-valued.

( )               + = ∑ - = k n N n x k X N n 2 1 cos ) ( ˆ1 0 π (4)
The envelope peaks, whose number and width are determined by the model order, will now be their frequency domain counterparts, the rectified waveform peaks. Thus, the model order has to be adjusted with respect to the temporal structure of the signal, and not to the formant structure of the spectrum.

TRUE ENVELOPE ESTIMATION

The true envelope estimator [START_REF] Villavicencio | On Cepstral and All-Pole Based Spectral Envelope Modeling with Unknown Model Order[END_REF] has been shown to outperform linear prediction [START_REF] Makhoul | Linear prediction: A Tutorial Review[END_REF] or cepstral methods such as discrete cepstrum [START_REF] Rodet | An Improved Cepstral Method for Deconvolution of Source-Filter Systems with Discrete Spectra: Application to Musical Sound Signals[END_REF] both in terms of accuracy and ease of model order selection. Recently the iterative procedure has been significantly improved such that the computational costs are in the similar to the costs of the Levinson recursion such that real time processing can be achieved [START_REF] Villavicencio | On Cepstral and All-Pole Based Spectral Envelope Modeling with Unknown Model Order[END_REF]. True envelope estimation is based on cepstral smoothing of the log amplitude spectrum and the resulting estimation can be interpreted as the best band limited interpolation of the major spectral peaks in such a way that the peak matching is maximized and inter-peak valleys are avoided.

Cepstral Smoothing

The real cepstrum is usually defined as the inverse Fourier transform of the log magnitude spectrum [START_REF] Villavicencio | On Cepstral and All-Pole Based Spectral Envelope Modeling with Unknown Model Order[END_REF], as shown in equation ( 5)

( ) ( )       = ∑ - = N kn j k X n x K k π 2 exp log ˆ1 0 (5).
Regarding the log magnitude spectrum as a signal, we can interpret each cepstral coefficient as a measure of the energy present in discrete frequency bands of that signal. Low-pass filtering the cepstrum (also called liftering) would result in a smoother version of the log magnitude spectrum, given by equation ( 6)

( ) ( )       - = ∑ - = N kn j k x w k C N n n π 2 exp 1 0 (6)
where C(k) is the smoothed spectrum (corresponding to the spectral envelope estimation) and w n is a low-pass window in the cepstral domain usually defined as [START_REF] Drullman | Temporal Envelope and Fine Structure Cues for Speech Intelligibility[END_REF] where n c is the cutoff quefrency. If we only want to represent the spectral envelope, discarding information about the partials we should set the cutoff quefrency below the period of the signal. One major drawback of this operation is that we discard spectral energy when setting cepstral coefficients to zero. The result is a smooth curve that is always below the peaks of the log magnitude spectrum. The true envelope estimator uses cepstral smoothing in an iterative procedure that aims at connecting the peaks of the log magnitude spectrum as explained below.

True Envelope

Let X(k) be the K-point DFT of the signal frame x(n) and C i (k) the smoothed spectrum at iteration i. The algorithm then iteratively updates the resulting spectral envelope A i (k) with the maximum of the original spectrum and the current spectral envelope

C i-1 (k) ( ) ( ) ( ) ( ) k C k X k A i i 1 , log max - = (8)
and applies cepstral smoothing to A i (k) to obtain C i (k). The procedure is initialized setting A 0 (k) = log|X(k)| and starting the cepstral smoothing to obtain C 0 (k).

Optimal Order Selection

The order of the cepstral representation of the spectral envelope is the number of cepstral coefficients we keep in the cepstral smoothing procedure, and as such is proportional to the fundamental frequency of the original signal. The optimal order should give a spectral envelope that follows the overall shape of the filter without representing the harmonic structure of the spectrum. In order to estimate the optimal order, we use the sourcefilter model and think of the spectrum as the result of the interaction of two components, represented by the source, an input signal that contains information about the frequencies of the partials and the filter that shapes the source spectrum. According to this model, the spectral envelope represents the filter that has been excited by the source. For near harmonic sources, the resulting spectrum will be quasi-harmonic. In terms of the interaction between source and filter, we can think of the resulting harmonics sampling the filter with a sampling rate that depends on the fundamental frequency of the source spectrum. According to the sampling theorem, we must sample the filter with at least twice the maximum frequency present in that signal. If we assume that the spectral envelope should not contain information about the harmonic structure of the spectrum, the maximum frequency present in that signal is the fundamental frequency F 0 , such that the related Nyquist frequency (assuming a sampling rate of Fs) is Fs = (2F 0 ). This formula provides a simple way of selecting the cepstral order because higher sampling frequencies would reveal (maybe partially) information about the harmonic structure of the spectrum and lower sampling frequencies would smooth out the spectral envelope, not revealing information about the (formant) peaks. We can therefore postulate the near optimal order of the cestrum given only that the maximum frequency difference between two spectral peaks that carry envelope information is known. If the difference between those peaks is Δ F then the cepstral order should be 5 . 0 , 2

ˆ0 = = ∆ = α α F F F O S F S (9).
While the optimal order, that is the order that provides an envelope estimate with minimum error, depends on the specific properties of the envelope spectrum, the order selection according to equation ( 9) is reasonable for a wide range of situations and the resulting error is generally rather close to the one obtained with the optimal order.

TRUE AMPLITUDE ENVELOPE (TAE)

Ideally, the amplitude envelope should be a curve that outlines the waveform, following its general shape without representing information about the harmonic structure. One of the most challenging aspects of this problem is that we are looking for a curve that is smooth during rather stable regions of the waveform, while being able to react to sudden changes (such as percussive onsets) when they occur. Here, we propose to use a dual of true envelope in the time domain. The time domain signal is subjected to the algorithm instead of the Fourier spectrum. In this way, the amplitude envelope is expected to match the amplitude peaks corresponding to the period of the waveform more closely than the previously introduced methods. The idea behind TAE is to mimic the structure of the spectrum with the time-domain signal to be able to apply the true envelope method directly. The basic steps to estimate the TAE are as follows. First we obtain a rectified version of the waveform (so that that are no negative amplitudes), next we zero-pad the rectified waveform to nearest power of two (thus mimicking the DFT), and then we finally add a time-reversed version of the zero-padded rectified waveform to represent the negative frequencies. Before estimation, we still need to exponentiate the amplitudes because true envelope supposes that we are fitting a smooth curve to the log magnitude spectrum. The result is illustrated in Fig. 1. The last step is the application of the true envelope estimation technique to obtain the true amplitude envelope (TAE), represented as a solid line outlining the rectified waveform.

It is important to notice that the peaks of the waveform do not carry the same information as the spectral peaks. Each peak of the spectrum corresponds to a partial, such that for quasi-harmonic spectra the separation between spectral peaks is F 0 . On the other hand, in only one period, the peaks of the half-wave rectified waveform generally contain information about all the frequencies contained in that signal (depending on their phases). Therefore, the time-domain counterpart of the near optimal order selection must take into account only the period of the waveform, instead of all rectified peaks. The optimal order is now directly proportional to the fundamental frequency of the waveform, instead of inversely proportional when using true envelope in the spectral domain because the separation of the spectral peaks Δ F is now represented by Δ T and given by the period of the signal T 0 for a half-wave rectified waveform (hwr) as equation ( 10 A full-wave rectified (fwr) version would present twice as many main peaks, requiring half T 0 , or α=1.

EVALUATION

This section compares the amplitude envelopes obtained with LPF, RMS, FDLP and TAE for a sustained and a percussive musical instrument sounds and a speech utterance. Figure 2 shows a halfwave rectified (hwr) version of a bass clarinet, a clavinet and speech utterance waveforms and the amplitude envelope estimates. We are looking for the estimation that best fits the model waveforms, following the amplitude evolution by matching the main peaks while avoiding ripples in more stable parts. It should be noted that FDLP was normalized and scaled to the maximum of the hwr waveform and that RMS was also low-pass filtered to eliminate the ripples during mostly the more stable parts. Upon close inspection, Figure 2 shows that TAE renders the best fit, closely following the peaks without ripples. TAE is also very responsive to sudden changes in amplitude, as can be seen for the clavinet sound. We should notice that the fit depends largely on the order. For isolated notes, the optimal order is expected to be the same throughout, but for speech the situation is not so simple because the F 0 changes dynamically. For speech utterances we can use the mean F 0 (supposing it does not change a great deal).

CONCLUSIONS AND FUTURE PERSPECTIVES

We proposed the true amplitude envelope (TAE) estimation technique as an improvement to the classical methods found in the literature. We have shown that TAE outperforms them by providing a smooth function that approximately matches the main peaks of the waveform avoiding under or over estimation due to the near optimal order selection based on the period of the signal. One important feature of TAE due to its cepstral nature is the ability to represent well sudden changes in the waveform while avoiding ripples during more stable parts. Future perspectives of this work could include evaluating the perceptual impact of the resulting envelopes with listening tests, developing an objective measure to evaluate the results and order estimation for signals with time varying F 0 .

Fig. 1 .

 1 Fig. 1. True amplitude envelope estimation. The figure shows the half-wave rectified and zero padded (zp) version of the waveform with its time-reversed counterpart used in the true amplitude envelope estimation method.Time

ACKNOWLEDGEMENTS

This work is supported by the Brazilian Governmental Research Agency CAPES (process 4082-05-2).


Marcelo Freitas Caetano

Xavier Rodet

SOUND MORPHING BY FEATURE INTERPOLATION

Keywords: Sound morphing, sonic features, musical instrument sound, timbre space

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Sound morphing has been used in music compositions [START_REF] Wishart | On Sonic Art[END_REF], [START_REF] Mcnabb | Dreamsong: The Composition[END_REF], [START_REF] Harvey | Mortuos Plango, Vivos Voco: A Realization at IRCAM[END_REF], in synthesizers [START_REF] Haken | Timbre Morphing of Sounds with Unequal Numbers of Features[END_REF], and even in psychoacoustic experiments, notably to study timbre spaces [START_REF] Moorer | Perceptual Evaluations of Synthesized Musical Instrument Tones[END_REF]. When morphing musical instrument sounds, we usually want to obtain hybrid sounds that are perceptually intermediate across timbre dimensions, such that the intermediate sounds would correspond to hybrid instruments between source and target. A challenging aspect is to control the transformation with a single parameter α, called morphing or interpolation factor [START_REF] Rodet | Automatic Timbral Morphing of Musical Instrument Sounds by High-Level Descriptors[END_REF], as illustrated in Fig. 1 for images. Figure 1 shows that, ideally, we want the morphing factor α to control perceptually related features of the transformation, such that the morph should be perceptually halfway when α = 0.5, for instance. Most morphing techniques proposed in the literature use the interpolation principle, which consists in interpolating the parameters of the model used to represent the sounds regardless of features [START_REF] Haken | Timbre Morphing of Sounds with Unequal Numbers of Features[END_REF], [START_REF] Covell | Automatic Audio Morphing[END_REF], [START_REF] Haken | Sinusoidal Modeling and Manipulation Using Lemur[END_REF]. In this work, parameter refers to coefficients from which we can resynthesize sounds, while feature refers to coefficients used to describe or identify a particular aspect of a sound. Usually, we cannot resynthesize sounds directly from feature values. The basic idea behind the interpolation principle depicted in Fig. 2 is that if we can represent different sounds by simply adjusting the parameters of a model, we should obtain a somewhat smooth transition between two (or more) sounds by interpolating between these parameters These authors often conclude that the linear interpolation of the parameters of known models does not correspond to linearly varying perceptually relevant features [START_REF] Covell | Automatic Audio Morphing[END_REF]. If we want the result to sound perceptually intermediate, we need to develop techniques to interpolate perceptually motivated features. For such, we adopted the morphing by feature interpolation principle, depicted in Fig. 2. Ideally, we would like to be able to interpolate in the feature space and retrieve the set of parameters that correspond to the interpolated feature values. Unfortunately, this is a notoriously difficult problem to solve since most features commonly used do not allow direct inversion for resynthesis, particularly when the features are correlated to perceptual characteristics of sounds [START_REF] Mcadams | Acoustic Correlates of Timbre Space Dimensions: A Confirmatory Study Using Synthetic Tones[END_REF], [START_REF] Mcadams | Caractérisation du Timbre des sons Complexes. II: Analyses Acoustiques et Quantification Psychophysique[END_REF], [START_REF] Winsberg | Perceptual Scaling of Synthesized Musical Timbres: Common Dimensions, Specifities and Latent Subject Classes[END_REF]. In this work, we describe techniques to obtain morphed sounds whose values of features are as close as possible to the interpolated feature values. Ideally, we want the feature values to vary linearly when the morphing factor varies linearly. Most features we use (log attack time, temporal centroid, spectral centroid, spread, skewness, kurtosis [START_REF] Peeters | A large set of audio features for sound description (similarity and classification) in the CUIDADO project[END_REF]) are acoustic correlates of timbre dimensions obtained by perceptual studies, such that sounds whose feature values are intermediate between two would be placed between them in the underlying timbre space used as guide. We measure the perceptual impact of the morphed sounds directly by the feature values, using them as an objective measure with which to evaluate the results. Thus we consider that the morphed sounds change perceptually linearly when the corresponding feature values vary linearly.

The next section introduces timbre spaces obtained in psychoacoustic experiments and the most salient dimensions of timbre perception unveiled. Then, we present the sound model we developed, which includes amplitude and spectral envelope estimation, temporal segmentation, and sinusoidal and residual modeling. Finally, we explain the feature interpolation techniques that we developed specifically for musical instrument sounds.

ACOUSTIC CORRELATES OF TIMBRE SPACES

Multi Dimensional Scaling (MDS) techniques figure among the most prominent when trying to quantitatively describe timbre.

McAdams [START_REF] Winsberg | Perceptual Scaling of Synthesized Musical Timbres: Common Dimensions, Specifities and Latent Subject Classes[END_REF] and Handel [START_REF] Handel | Timbre perception and auditory object identification[END_REF] independently propose comprehensive reviews of the early timbre space studies. Grey [START_REF] Moorer | Perceptual Evaluations of Synthesized Musical Instrument Tones[END_REF] investigated the multidimensional nature of the perception of musical instrument timbre, constructed a three-dimensional timbre space, and proposed acoustic correlates for each dimension. He concluded that the first dimension corresponded to spectral energy distribution (spectral centroid), the second and third dimensions were related to the temporal variation of the notes (onset synchronicity). Krumhansl [START_REF] Krumhansl | Why is Musical Timbre So Hard to Understand?[END_REF] conducted a similar study using synthesized sounds and also found three dimensions related to attack, synchronicity and brightness. Krimphoff [START_REF] Mcadams | Caractérisation du Timbre des sons Complexes. II: Analyses Acoustiques et Quantification Psychophysique[END_REF] studied acoustic correlates of timbre dimensions and concluded that brightness is correlated with the spectral centroid and rapidity of attack with rise time in a logarithmic scale. McAdams [START_REF] Winsberg | Perceptual Scaling of Synthesized Musical Timbres: Common Dimensions, Specifities and Latent Subject Classes[END_REF] conducted similar experiments with synthesized musical instrument timbres and concluded that the most salient dimensions were log rise time, spectral centroid and degree of spectral variation. More recently, Caclin [START_REF] Mcadams | Acoustic Correlates of Timbre Space Dimensions: A Confirmatory Study Using Synthetic Tones[END_REF] studied the perceptual relevance of a number of acoustic correlates of timbre-space dimensions with MDS techniques and concluded that listeners use attack time, spectral centroid and spectrum fine structure in dissimilarity rating experiments. Here we should notice that most MDS techniques suppose that the underlying space is orthogonal and metric, which means that the dimensions are independent and the notion of distance is defined. Therefore, shifting linearly from a source to a target sound in such space would correspond to a perceptually linear change across all dimensions and would also result in linear variation of the values of the correlates of each dimension.

The Features Used as Guides

We included temporal and spectral features related to the acoustic correlates of timbre spaces presented earlier [START_REF] Peeters | A large set of audio features for sound description (similarity and classification) in the CUIDADO project[END_REF] and are supposed to capture the most perceptually salient dimensions of the timbre spaces proposed, namely, the attack time and the distribution of spectral energy. The temporal features we use are the log attack time and the temporal centroid. The spectral shape features we take into consideration are spectral centroid, spread, skewness and kurtosis.

MODELING PERCEPTUALLY RELEVANT FEATURES

The morphing technique we developed is depicted in Figure 4, and consists of two basic steps, temporal segmentation and consequent temporal alignment followed by spectral modeling and spectral morphing. The temporal processing steps are represented by light grey blocks and spectral processing by blocks with white background. In general terms, the sound model consists in a temporal segmentation step, amplitude envelope estimation, sinusoidal plus residual modeling and spectral envelope extraction.

Temporal Modeling

Temporal modeling comprises the temporal segmentation step and amplitude envelope estimation. The temporal segmentation, shown in Fig. 5, consists in automatically estimating the boundaries of four perceptually important regions, namely, the attack, the transition, the steady-state or sustain, and the release using an automatic segmentation technique we proposed elsewhere [START_REF] Burred | Automatic Segmentation of the Temporal Evolution of Isolated Acoustic Musical Instrument Sounds Using Spectro-Temporal Cues[END_REF]. From these estimations, we calculate the length of each region. The most salient descriptor here is log attack time. Amplitude envelope estimation is performed with the true amplitude envelope technique we developed [START_REF] Burred | Automatic Segmentation of the Temporal Evolution of Isolated Acoustic Musical Instrument Sounds Using Spectro-Temporal Cues[END_REF], based on cepstral smoothing. The interpolated amplitude envelope modulates the spectral frames of the morphed sound, so we use the temporal centroid as its feature.

Spectral Modeling

The basic spectral modeling technique used is harmonic sinusoidal modeling plus noise residual [START_REF] Serra | Musical Sound Modeling with Sinusoids Plus Noise, Musical Signal Processing[END_REF]. In order to independently manipulate the spectral shape, we model the amplitude of the partials using a spectral envelope model instead of the amplitudes output by the sinusoidal analysis. Fig. 6 illustrates both representations of the sinusoidal model. In part a) we show the traditional sinusoidal representation, where each partial has a frequency and associated amplitude tied together. In part b), however, the spectral envelope controls the global spectral shape independently from the frequencies of the partials. For every frame of both the sinusoidal and noise residual we calculate the spectral envelope using true envelope [START_REF] Rodet | Efficient Spectral Envelope Estimation and its Application to Pitch Shifting and Envelope Preservation[END_REF]. The sinusoidal component is modeled as shown in part b) of Fig. 6, and the noise residual is modeled as white noise filtered with the spectral envelope estimated from each frame of the noise residual.

MORPHING STRATEGY AND EVALUATION WITH FEATURE VALUES

The morphing steps are as follows: temporal alignment, spectral envelope morphing and amplitude envelope morphing. We use the values of features as a guide to measure the perceptual impact of each step. Ideally, we want the features to vary linearly when the morphing factor varies linearly.

Temporal Alignment

The first important step in morphing is the temporal alignment of perceptually different regions, such as the attack, characterized by fast transients, and the sustain part, much more stable. We cannot expect to attain good results if we combine a sound that has a long attack with another sound with a short one regardless of the differences. The region where attack transients are combined with more stable partials will not sound natural. So, with this in mind, to achieve a more perceptually seamless morph, we need to temporally align these regions so that their boundaries coincide.

For each sound, we measure the length of each region (labeled with letters) by computing the time difference using the markers (numbers). The length of the attack is represented in logarithmic scale. For the other regions the representation is linear. Then, we interpolate between the lengths of the regions according to equation ( 1) to obtain their corresponding lengths in the morphed sound. The interpolated lengths are represented by a letter that stands for the region and subscripts indicating both sounds, e.g. S 12 for the sustain as shown in equation ( 1)

[ ] 2 1 12 1 S S S α α - + = (1)
where S 1 represents the length of the sustain of the first sound and S 2 of the second. The stretch/compress factors R S1 for the first sound and R S2 for the second are calculated as in equation ( 2) Finally, we simply time stretch/compress each region by the corresponding ratio, for instance, S 1 by R S1 , etc. Since we interpolate the attack time logarithmically, the attack times (perceived on a logarithmic scale [START_REF] Mcadams | Acoustic Correlates of Timbre Space Dimensions: A Confirmatory Study Using Synthetic Tones[END_REF], [START_REF] Winsberg | Perceptual Scaling of Synthesized Musical Timbres: Common Dimensions, Specifities and Latent Subject Classes[END_REF]) will be perceived as varying linearly under the time-alignment transformation. The other regions will also be properly aligned and ready to be morphed in the spectral domain.

Spectral Envelope Morphing

Following the morphing by feature interpolation principle, the objective of the spectral envelope morphing step is to obtain a morphed spectral envelope that has intermediate formant peaks and intermediate values of spectral shape features. In order to achieve this, the ideal would be to be able to interpolate the spectral feature values and invert this representation to obtain spectral envelope parameters corresponding to the interpolated feature values. One difficulty this approach poses is that there is no known analytic inversion from the chosen features. Instead, we will study which spectral envelope representation leads to linearly varying values of spectral shape features when its parameters are linearly interpolated. The representations investigated are the envelope curve (ENV) [START_REF] Haken | Timbre Morphing of Sounds with Unequal Numbers of Features[END_REF], [START_REF] Haken | Sinusoidal Modeling and Manipulation Using Lemur[END_REF], line spectral frequencies (LSF) [START_REF] Paliwal | Interpolation Properties of Linear Prediction Parametric Representations[END_REF], cepstral coefficients (CC) [START_REF] Covell | Automatic Audio Morphing[END_REF] and dynamic frequency warping (DFW) [START_REF] Meyers | Morphing Spectral Envelopes using Audio Flow[END_REF]. These correspond to the main sound morphing methods proposed in the literature.

Figure 7 shows the source and target envelopes in solid lines and nine intermediate envelopes corresponding to linearly varying the interpolation factor α by 0.1 steps in dashed and dotted lines; on the right, we see the associated values of the spectral shape descriptors for each step. Figure 7 confirms that for this case interpolating envelope curves does not account for formant shifting and most spectral shape descriptors do not vary in a straight line. Figure 7 also shows that the linear interpolation of cepstral based envelope representations like Slaney [START_REF] Covell | Automatic Audio Morphing[END_REF] proposes neither shifts the formants nor results in linear variation of descriptors. The same applies for the DFW based spectral envelope morphing proposed by Ezzat [START_REF] Meyers | Morphing Spectral Envelopes using Audio Flow[END_REF]. On the other hand, LSFs behave fairly well under both constraints in this case just like Paliwal [START_REF] Paliwal | Interpolation Properties of Linear Prediction Parametric Representations[END_REF] states for LSFs. In conclusion, in our model we adopt LSFs as the most suitable parameters to represent and interpolate the spectral envelopes.

Interpolation of Partial Frequencies

Since we are morphing quasi-harmonic musical instrument sounds, there exists a direct one to one correspondence between the partials of both sounds. All we have to do is make sure we will be able to fuse them into a single percept without losing the inherent spectral variations. Here we propose to morph quasi-harmonic musical instrument sounds with the same pitch, so that the partials have a one to one correspondence and no pitch shift is required. Since the spectral shape and form are morphed separately with the spectral envelope, we simply interpolate the partials frequency values to account for frequency fluctuations (jitter, shimmer), inharmonicity and other temporal features that are encoded in the frequency variation with time.

Morphing the Amplitude Envelope

Morphing the amplitude envelope is similar to the spectral envelope because the techniques we developed for estimating the amplitude envelope are inspired by spectral envelope estimation techniques [START_REF] Burred | Automatic Segmentation of the Temporal Evolution of Isolated Acoustic Musical Instrument Sounds Using Spectro-Temporal Cues[END_REF]. Also, the temporal centroid is the time-domain analogous of the spectral centroid [START_REF] Peeters | A large set of audio features for sound description (similarity and classification) in the CUIDADO project[END_REF], and as such, its values behave in the same fashion under the same transformations. Fig. 8 confirms that the temporal centroid varies fairly linearly even when we interpolate the amplitude envelope curves directly.

Morphing the Noise Residual

We morph the spectral envelopes of the residual noise signal and synthesize a morphed residual by filtering white noise with it and mixing it into the morphed sinusoidal component. The noise residual is obtained simply by filtering white noise with a timevarying filter obtained from the interpolation of the spectral envelopes extracted from the residual of the original sounds. We interpolate the LSFs used to represent the spectral envelope of the noise residual just like for the sinusoidal part.

CONCLUSIONS AND FUTURE PERSPECTIVES

In this work, we describe techniques to morph salient timbral dimensions of quasi-harmonic musical instrument sounds using features as guides. The features we use are acoustic correlates of timbre dimensions obtained in psychoacoustic studies, such that sounds whose feature values are intermediate between two would be placed between them in the underlying timbre space. So, interpolating the feature values becomes the goal itself to render the results more perceptually linear. For such, we introduced the concept of morphing by feature interpolation and developed a model and techniques to apply it. The sound model aims at allowing independent manipulation of perceptually meaningful features of sounds related to salient timbre dimensions, such as spectral and amplitude envelopes, attack time, among others. We measured the perceptual impact of the morphed sounds directly by the feature values. We consider that the morphed sounds change perceptually linearly when the corresponding feature values vary linearly. Future perspectives could include vibrato and tremolo modeling and treatment, and extending the technique to inharmonic sounds, which would probably require a technique to find correspondences between the partials of both sounds. Sound examples on http://recherche.ircam.fr/anasyn/caetano/icassp2011.html.

Fig. 2 .

 2 Fig. 2. Morphing using the interpolation principle.

Fig. 4 .

 4 Fig. 4. Depiction of the general steps in our morphing procedure. The blocks represent temporal and spectral feature extraction and processing.

Fig. 5 .

 5 Fig. 5. Amplitude/Centroid Trajectory model. The figure shows the regions indicated by letters and their boundaries indicated by numbers.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Spectral representation of partials. The figure shows the traditional sinusoidal representation with the frequency values and amplitudes tied to each other in part a). Part b) depicts our representation, where the amplitudes of the partials are represented independently with a spectral envelope model.b) a)Amplitude

Fig. 8 .

 8 Fig. 8. Amplitude envelope morphing guided by high-level temporal shape descriptors. The figure shows the variation of the temporal centroid value when interpolating between the parameters of the amplitude envelope. The curves are shown on the left and the corresponding descriptor variation on the right.

Fig. 8 .

 8 Fig. 8. Amplitude envelope morphing guided by high-level temporal shape descriptors. The figure shows the variation of the temporal centroid value when interpolating between the parameters of the amplitude envelope. The curves are shown on the left and the corresponding descriptor variation on the right.

ACKNOWLEDGEMENTS

This work is supported by the Brazilian Governmental Research Agency CAPES (process 4082-05-2).