Biofilm early stages of growth and accumulation theoretical model

A theoretical model to translate the evolution over time, in early stages, of growth and accumulation of biofilm bacterial mass is introduced. The model implies the solution of a system of differential-difference master equations. The application of an algorithm like Miller ́s tree term recurrence, already known for Bessel functions of first kind, allows an exact calculation of the solutions of such equations, for a wide range of parameters values and time. For biofilm model a five term recurrence is deduced and applied in a backwards computation. A suitable normalisation condition completes the reach of the solution.


Introduction
Discrete mathematical models in space, and continuous in time, arise in many different contexts of Physics, Chemistry and Engineering. Ever since, consideration of more sophisticated exact mathematical solutions, for transient regimes, invariantly leads, in known bibliography, to numerical and approximated methods application. Such approach neglects all the theoretical information contained in an exact solution, which usually consists of special functions or, more frequently, complex combinations between them (in particular an infinite sum of terms with Bessel functions of first kind, modified or not).
As a matter of fact the accurate calculation of special functions seems far from being an inviting task in several fields of exact sciences, namely Operational Research [1,2], Teletraffic Theory [3] and Manufacturing Systems [4]. The commentary in [2] and subsequent ones, with the same inspiration, many years after [5], and even nearly nowadays [6] corroborate such difficulties. Functions similar to that of transient M/M/1 queue also came into play in Physical-Chemistry (nucleation, polymerization, crystallization, ... , aggregation processes in general), in Radar Theory [7], in Computer Networks Theory [8], and Congestion phenomenon [9].
From those exact solutions little information is brought to light, with the exception of long time behaviour (steady state) which is a limit more easy to deduce. If transient regime is considered, having already an exact solution, a formalism of type Fokker-Planck is usually chosen to approximate such regime. The purpose is to mimic a process discrete in space and continuous in time by a process continuous both in space and time.
Doing things like that the initial infinite system of differential-difference equations is transformed in an equation with partial derivatives (Fokker-Planck equation). On the other hand special functions, that are the most important part of the solution, are eliminated from all the analysis, and so we lose all their formal beauty.
Formalism conducting to Fokker-Planck equation started when Frenkel brought Becker-Doring theoretical work on nucleation kinetics to "its modern form", as Goodrich said [10]. Seems that such "modern form", or others similar approximations, has been generally accepted in many domains of scientific knowledge already cited.
In recent years a theoretical model to account quantitatively the evolution over time of growth and accumulation of biofilm biomass lead to a few similar exact solutions [11,chs. IV and VI]. In particular, the simplest version of all variants analysed there, lead to an analytical solution at least so complex as that of the M/M/1 queue [11, pp. 67 and 70].
To compute the solution with accuracy one must deduce suitable homogeneous recurrence relations of a part of it, containing, at least, the infinite sum of terms with the Bessel functions. This has been in fact the case for the aforesaid simplest version of the biofilm model. A suitable adaptation of the already known three terms recurrence (TTR) Miller algorithm [12] for the functions ( ) and ( ) to this new context in which recurrences have now more than tree terms has been done.
Miller applied backward recurrence "for evaluating sequences of functions { fk } when the recurrence connecting successive members was unstable for increasing k" as published Thacher [13]. Applying an algorithm like Miller´s (TTR) [12], we performed an exact calculation and graphic representation of such kind of solutions, with a five term recurrence, also applied in backwards calculation. Consequently we have applyed, in the simplest version of our biofilm model, a generalization of that algorithm including the whole infinite sum in functions { fk } and not only one Bessel function. So the recurrences obtained have more than tree terms, as expected according to the additional complexity of the functions.
The early stages of growth and accumulation of biofilm included in that model, was the subject of [11] jointly with the search of how to deal when facing theoretical models near the border of mathematic intractability.
Accordingly, those early stages can be modeled by several more complicated versions than the one deduced in this work, which is only a first approach.
Besides this first version wich we call "Mono-layered concentrated growth" kinetics, we, also, for sake of completeness, described in [11] how to reach solutions for three more approaches: Zero order kinetics, "Blackman inspired" kinetics and "Averaged" first order kinetics, all known in biotechnological literature.
Those approaches are, however, only particular cases of the main task, which is the mathematical deduction of the exact solution for Monod kinetic law. We have already reach such general version of the model [11, Chap. VI].
As conclusion a Bessel functions framework seems to be appropriate for modeling biofilm formation and growth dynamics. The ideas and the model described here are only the beginning of the modeling possibilities offered by that mathematic framework.
In Queuing Theory utilization of State-Transition-Rate Diagrams is a customary tool for description of a particular queue dynamics under analysis. We also constructed the equivalent diagram for biofilm model, lightning this way an "isomorphism" between the two contexts [11, p. 182]. Accordingly with this circumstance, a careful analysis of Clarke classical solution [14] and our previous work [11] allows to put in equation and reach the transient solution of many more queues [11,Chap. VII,. Such unifying way will allow, with a unique methodology, obtain exactly solutions in many models. Transient phases in Queuing Theory is a subject that we will leave for future work.

Skope and structure of this work
The biofilm theoretical model is briefly described in Part I. The mathematical resolution is accomplished in Part II, the deduction of an algorithm similar to Miller´s three term recurrence is deduced in Part III and some numerical simulations are presented in Part IV. Conclusions and future work are in Part V.

Theoretical model brief description
A brief description of the model can be made with figures 1, 2 and 3. The model postulates a layered structure to describe that the thickness of the biofilm is not uniform as in Fig. 1.  It is also considered a gradient of nutrient concentration, from the outermost layer into the biofilm, along L layers. The nutrient concentration profile is postulated to decreases linearly, from layer to layer, from the outside to the inside, until it reaches zero after (L + 1) layers. This linearization is only a first approach in the context of this most simple version of the model. In Fig. 3, it is L = 4. active (colored from red to yellow) and inactive (in white).
In active population red means the highest concentration of substrate, yellow the lowest one and in inactive population white means that there is no substrate In areas where > , all the layers below the first L from above are inactive and do not consume nutrients. The bacteria in the layers in which exists nutrient are active and can reproduce, thus contributing to the increase in biomass. The nutrient concentration in each layer determines the specific growth rates of bacteria in that layer.
The partial processes that contribute to the formation and growth of biofilm are: 3 -The specific growth rates of bacteria in the active layers, quantified by = á . · ( ) in which á . is the maximum specific growth rate and is the remaining factor that depends on the nutrient concentration at the considered layer.

-
The maximum specific growth rates are decomposed in the directions parallel and perpendicular to the solid support thus getting defined two specific growth rates: → = µ á . → · ( ) and ↑ = µ á .

↑ · ( )
More minor details can be found in [11,Chaps. I,II,III]. Accounting for all these partial processes we obtain a system of differential-difference master equations whose solution will be the functions = ( ) for ≥ 0.

Part II
Theoretical model resolution

Governing equation system definition
We start the resolution of our model writing the total system of differential-difference master equations obtained in [11, Chap. III, Box 5].
The system reads, from (1) to (4): -number, from 1 to , and from top outermost layer to down innermost still active layer in any area [11, Chap. III, pp. 25].
All terms with common S are colected together and are disposed in increased order ( ) in the right hand sides. We observe that the system structure consists in, 6 of 43 : one differential equation The dimensionless variables and parameters that we choose are defined from (5) to (11): Before insert the dimensionless parameters it is convenient put in evidence the kinetic factor whenever it occurs to facilitate the substitution of . This is done from (12) to (15).
As consequence of (5) and (6) the derivatives on the left hand sides assume the form: The next step is to divide all the equations by µ á . → + µ á .
↑ · · and we get immediately the dimensionless system, equations (17) to (20). Besides, initial condition (21) is part of the system, consequently also must be included.
The initial condition must, of course, be: (0) = 1, and (0) = 0 ˃0 (21) ( at the beginning (τ=0) al the solid support is bare) 5. General model resolution strategy: concept of "Mathematical meeting point" This system of equations is composed by the so-called differential-difference equations. It is a traditionally difficult domain and suitable theory for this subject can be found in Pinney [15], Bellman and Cooke [16] and El´sgol´ts and Norkin [17] and certainly other more recent works. Also formally similar problems arise in Operational Research literature, in particular when transient states are analysed in Queuing Theory. In the present case the system has an analytical solution. However that solution can´t be represented, for the most general version, by a closed form, like a final expression = ( ), and we will only solve the model for a simplified case and extract observations from that resolution. Such observations will later allow to establish the method for more general cases. A formalism similar to that of Clarke ([14], [18]) will be our working tool. Suitable modifications will be introduced to adapt such solving way to the particular details of our case. Clarke formalism can also be found described in Saaty [19] and Srivastava and Kashyap [20]. We will use the following strategy, in two steps, being only the first one the scope of this work. The second one will be next future work: First step: We must solve the system for = 1 and = = ⋯ = = ⋯ = . That is, for the case where there is only one active layer and erosion is of equal intensity in all the surface fractions ( ≥ 1) . This one will be named the simplified version of the model. More specifically we can name this version the "Mono layered concentrated" growth kinetics variant, already defined in [11, Chap. III]. Second step: analyse the solution obtained in the first step, inspect their mathematical structure, and construct a suitable Guess Solution that will reduce all the solution of the complete system to an algebraic problem. The algebraic problem defined in this second step falls in the context of algebraic difference equations and is somewhat cumbersome. However it will be demonstrated that the complete resolution is possible. The aforesaid Guess Solution contains an infinite series of Bessel functions, each one multiplied by a respective labelled factor. All those factors must then be analytically determined on a recurring suitable scheme. The system of the simplified "Mono layered concentrated" growth kinetics version, corresponding to the first step, is easily obtained from (17) to (21) inserting the conditions: = = ⋯ = = ⋯ = (22) ... and, = 0 , … … … ˃ 1 … … … ≥ 2 ( 23) Condition (23) is equivalent to impose = 1 because all kinetic parameters are annulled for ≥ 2. We also observe that (23) implicates the sameness between (18) and (19). Consequently, the system for the simplified version is: We have already defined, as work context, [[11, p. 19, figs. 6 to 11] the five specific growth rate kinetics: -Monod kinetics -Zero order kinetics -"Blackman inspired" kinetics -"Averaged" first order kinetics -"Mono-layered concentrated growth" kinetics 9 of 43 When = 1 the first four (Monod, Zero order, "Blackman inspired" and "Averaged" first order) reduce in complexity and becomes equal to "Mono-layered concentrated growth" kinetics model. So we can conclude that to solve this last model is equivalent to solve the most simplified version of all the other four. We can consider "Mono-layered concentrated growth" kinetics as a sort of "mathematical meeting point" of all the kinetics under our analysis.
6. "Mono layered concentrated" growth kinetics: exact solution 6.1 Common initial steps: generating function formalism The first step to start the resolution is to make the substitution, This substitution aims to simplify the system by eliminating the term relative to index ( ) on the right hand side, in all the equations for ( ≥ 1). After simplification we get, from (29) to (31), With initial condition: Let now define the following generating function, ... similar to that one of Clarke ([14], [18]), ... but with additional coefficients which remain to be determined and with changed signs at and .
These modifications are justified because the system solved by Clarke is, Therefore this is a simpler system because and are null which is not the case of our system. This explains the need to include the coefficients . Another difference lies in the exchanged signs in the functions ( ) which explains the need to change the signs in and in the generating function. Like in Clarke ([14], [18]) here we also seek for a generating function, ( , ) such that obeys a differential equation like, Which is the so called Telegraph equation. The suitable deduction of coefficients will accomplish that goal.
First, ( , ) must be obtained from (33). We get, The order of partial derivations is indifferent, Now, taking the expressions for with ≥ 1 given by (30) and (31) we get, The usefulness of coefficients is now proved inspecting the second summation because we want to cancel it. The reason is that the exponents in ( − ) is two units lower than the corresponding order of functions and this summation cancelation would not be possible if all coefficients where equal to 1, like in Clarke´s generating function. Determination of coefficients is easy and follows straightaway.
Must be, We get immediately the solution, ... where, Now equation (29) and condition (32) will come into play. According to (32), at time ( = 0), the generating function must be ( , ) = (0) = 1. So, condition (32) now reads. in generating function formalism, With (49) the partial derivative can be easily obtained, Inserting here equation (29) one gets, ... and, What means, .. and, consequently, Equations (54) and (57) are now included in (53) and, in this way, the last condition for complete definition of the problem, in generating function context, is, Once deduced, the generating function ( , ) will conduct easily to all the functions ( ).
In fact from (49) one can conclude that, ... and, This last formula is easy to demonstrate, so we omit the details. Lastly, with (28), fractions ( ) are calculated from the functions ( ). Then we will reach the solution of the simplified version of the model which is also the "mathematical meeting point" for all the kinetic variants under analysis. A summary of what has been obtained so far in this generating function formalism, and will be used in the sequel, is, from (I) to (VII), ) ( , 0) = 1 (Initial time condition) (50) Telegraph equation (48) The function ( , ) is the Laplace transform of ( , ) and, Similarly, applying the transform to the right hand side, The right had sides in (62) and (64) must be equal, The initial time condition (50) after Laplace transform reads, ... which is equivalent to, The general solution of (65) is known from mathematical tables, The integration "constant" ( ) is only constant relatively to time and only depends on .
From (68), Substitution of (67) here leads to, Now is easy to particularize the general solution (68) for our case, inserting there this expression of ( ), An observation must be made about the change of by as integration variable. This has been done with the purpose of avoid confusion between that variable and the superior limit in the integral.

Trifurcation stage: continuous biofilm, patchy biofilm and border case
In the sequel, application of inverse Laplace transform leads back again to the generating function ( , ).
We must now consider tree hypotheses depending on the relative values of parameters and (1 − Φ). As bigger is , or as littler is (1 − Φ), more the process of erosion is important in their intensity, comparatively to the process of bacterial growth inside the biofilm. In mathematical context this is compulsory because if the sign in the exponentials is positive the inverse transform is one, if it is null is another and if it is negative is yet another one.
The functions in Cases 1 and 2 are respectively Bessel functions and Modified Bessel functions of first kind of order zero. Bessel function of first kind of order is defined by, ... and Modified Bessel function of first kind of order by, Here ( ) is the gamma function. When =1, 2, 3, ... then, In our case = 0, so ( + + 1) = ( + 1) = ! and we can write, ... and, As indicated before, in the first stages of biofilm formation erosion is seldom the most intense process. In most of situations growth as consequence of active bacterial reproduction is nearly the only relevant process.
This corresponds, in mathematical terms to the aforesaid Case 1. And we can classify the obtained biofilm as a "continuous biofilm".
However, patchy biofilms, also exists everywhere if erosion becomes strong enough in those initial times of development, leading to a steady state where bare areas of solid support exist in finite extension, even after a long time, and biomass remains only accumulated in colonies or clusters without a complete coverage of solid support area.
The designation "patchy biofilm" is opposed, or complementary to "continuous biofilm" in this physical classification. Consequently we observe that this mathematical formalism can translate this duality in biofilm structure, relating different biofilm morphologies with different mathematical functions. Such translation gets more visibility, principally when achieved formulas are graphically represented. Case 3 is the border case and the limit value of the bare fraction, ( ), is zero at infinite time , as expected.
Complete coverage of initial support bare area is rendered complete, after some time, in the case of con- Where: Solution for Case 2: patchy biofilm: Where: * = 2 · · ( + ) · ( + − 1) Solution for Case 3: border case between Cases 1 and 2: We provide the demonstration of these three solutions respectively in Appendix A, Appendix B and Appendix From these proof we concluded that to handle Case 1 solution (80), leads to a mathematical difficulty similar to handle Case 2 solution (81). And also leads to the same recurrent result.
Consequently we focus now our attention exclusively in continuous biofilm solution (80).
This solution is not easy to handle for numeric calculation and therefore also not easy for graphic and intelligible representation of the model. The solution obtained is complex because it consists of an infinite series of terms and the main difficulty lies on the fact that each term, in the infinite sum, contains a Bessel function of first kind, ( ).
It is not because the summation is infinite in the number of terms why it is difficult reach a rigorous calculation. The truth is that terms of bigger order will always diminish in absolute value, and this circumstance allows, sooner or later, to despise those upper order terms thus becoming possible to accomplish all the computation without error.
The main problem to solve refers to the rigorous calculation of each Bessel function not only because ex- This is an additional inconvenient question related to the lack of accuracy uniformity in data provided by known tables. Bessel functions of lower order (0, 1, 2, ... ) usually are tabulated with much more decimal significant places than the others (Abramowitz and Stegun, [21]). Thus, when one intends to calculate some more terms in the series, to avoid neglect still significant terms, coherence is lost because their accuracy is almost always inferior than that of the initial terms.
Consequently we must work only with computer tools paying attention to achieve the better exactitude disposable nowadays. We choose to compute in Excel spreadsheet not only because it stores numerical values as "Double Precision Floating Point", representing numbers accurate to something like 15 decimal places, but also because it allows to design the layout of the algorithms, over the spreadsheet with great personalized liberty. Another important advantage lies on the fact that one can "see" the contents of each cell and easily detect and correct errors and (or) improve the global efficiency of the entire algorithm.
Once the numerical values published in these tables are discarded, our attention turns to the theoretical method used, in past times, to construct such tables of Bessel functions.
The central and fundamental tool for such purpose was the Three Terms Recurrence (TTR) Miller´s algorithm (Miller,[22]).
The immediate obvious way to solve the problem, of handle solution (80) with adequate accuracy, is to apply TTR Miller´s algorithm, term by term, till reach a Bessel function order high enough for neglect the absolute value of such term and also the absolute values of all following terms.
Another, more efficient and elegant, alternative concerns with considering all the infinite series summation, and so all the function ( ), as a special function passible for application of an algorithm related to, but more general than, TTR Miller´s algorithm.
In more explicit words: functions ( ) will came into play relatively to this new algorithm in a similar way Bessel functions of first kind play their role relatively to the already known TTR Miller´s algorithm. Application of this new approach implies, as first step, the deduction of a corresponding homogeneous recurrence for the special functions ( ).
In the following section we will deduce that homogeneous recurrence and explain how can it be used for accurate computation of functions ( ), translating to this new method a similar reasoning to that one of the old TTR Miller´s algorithm method.
These two methods can be jointly, and independently, applied and so we will be able to validate mutually both.

Generalized Miller´s like algorithm construction.
The calculation of the model equation (80), will be done, in this section, applying a method alike that of Miller [22] but more general.
This method applies to the so-called Special Functions when one can construct recurrence formulas between them. Knowing the functions of lower orders we can calculate those of higher orders with a formula of progressive recurrence and vice versa with a regressive one.
However, if the Special Functions are decreasing when their order increases the application of the progressive recurrence will imply the use of the subtraction operation and the loss of significant numbers (S. Zhang and J. Jin, [23]) each time the recurrence computation applies, making the method numerically unstable.
This is the case for Bessel functions of first kind and therefore also the case for the functions ( ) at the initial values of ( ). But not only in the initial values of ( ) since there will always exist surface fractions ( ), in the largest orders ( ) that are in stage of increasing its value, having therefore a low value, since the kinetic process is an overall continuous accumulation of biomass.
Therefore, it is necessary to use a regressive formula in all values of dimensionless time ( ). That is the same to say that we have to compute the functions ( ) with the smallest orders ( ) depending of those ones with the biggest orders ( ).
Miller [22] solved the problem by assigning an arbitrary value to the function where the recurrence calculation starts, which is the highest order function (n) of all of them. By applying the homogeneous recurrence formula one gets the relative correct magnitudes between all the functions under calculation. After a suitable common normalization factor must be applied to all the functions of the sequence. The result of this multiplication is the exact value of all and each one of the functions.
The aforesaid suitable common normalization factor is defined by the known value, in the particular case under analysis, of an infinite series with form, The homogeneous recurrence formula to be constructed must, for the functions ( ), read like, ... where, ... and being ( ) obtained depending of ( ), ( ), … , … and ( ). We call (84) a ( + 1)-terms homogeneous recurrence.
Let´s deduce now the homogeneous recurrence, and for start, let define the auxiliary functions ( ), for ... where coefficients are defined like they have been defined in the model equation (80), Following a brief revision from Appendix A, coefficients are, in their turn, defined by equations ( 19(0)) and ( 29), And, also from Appendix A, consequently the coefficients obey equation ( 20), ... where, = · ( + − 1) Let use functions ( ) in the form, Now we properly multiply (89(1)) and (89 (2)) in such a way that, by associating the same Bessel function in the three equations (89(0)), (89 (1)) and (89 (2)), equality (A20) can be used thus eliminating the infinite terms summation: Adding the three equations (90(0)), (90(1)) and (90 (2)) we obtain, Observing that, from = · !, and also 19 (0) , 19 (1) , ( 22) and ( 23), the following three equalities are valid, ... we are now able to insert they in (91) , and therefore getting in it a great simplification, Let now add 2 in the discrete order ( ) at (95), ... and call into play the identities, ... which must be substituted in (95) and (96), The wanted regressive recurrence is now reached easily from (103), taking (104) in consideration, What logically follows must be to relate functions ( ) with functions ( ) and then insert that relation in The last two regressive recurrences, 108(0) and 108( ) can be clustered in only one, which reads like, General regressive recurrence for functions θ (τ) (biofilm early stages formation and growth model) Where: Remark: The designation "General" for this recurrence means that it is valid in the three cases analysed in the Appendixes A, B, and C, in Part II: Our goal is to apply correctly alike ideas to those applied by Miller [22] but now with this five terms recurrence (FTR) instead of the already known three term recurrence (TTR) classical algorithm. In the same formula, (108( )), a value of 1 must be assigned to ( ) if we follow strictly Miller´s criterion (Miller [22]). However we introduce in this stage an improvement consisting in assign to ( ) a small positive value, , but not necessarily equal to 1, and preferably less than 1. This allows a more con- The second recurrent computation consists in apply once again (108( )), but diminishing 1 to all orders ( ± ). We designate by (108( − 1)) such "new" equation, Now we insert the former values of ( ) = 0, ( ) = 0, ( ) = and ( ) = and obtain the non normalized value of ( ), coherently designated by .
This procedure is done repeatedly ( − 2) more times till the non normalized value of ( ), designed , is calculated.

Numerical examples
The model parameters are ≥ 1, > 0 , 0 ≤ ф ≤ 1 and ≥ 0 but, in this work, we limited the analysis to the simplified version with = 1.
In Fig. 4 the possible values of ф and δ divide the graph into two zones: the triangle ΔACI, corresponding to continuous biofilm, and the zone between ф = 0 and ф = 1 and above line CI, corresponds to patchy biofilm. Borderline case corresponds to points in line CI. With the indicated points, numerical simulations can be made in order to obtain a sufficiently enlightening output. For each point, it seems reasonable to vary ρ between 0.001 (near zero) and 10 (much above 1). We will illustrate the model, for continuous biofilm, with points A and E and, for patchy biofilm, with point D. In Fig. 5, 6 and 12, the functions = ( ) are represented for points A, E and D. In addition, Tables 2   and 3 contain numerical values that also illustrate the results obtained. Table 2 refers to the examples in Fig.   5 and 6, and Table 3 to those in Fig. 12.
In each one of the Fig. 5, 6 and 12 there are seven graphs corresponding, from left to right, and top to bottom, to the values of = 0.001, 0.01, 0.1, 0.2, 1, 5 and 10. Numerical values were calculated using the algorithm described in Part III and have 12 decimal places.
In Fig. 5 (ф = 0 and = 0) the only processes are the arrival and attachment flow and the growth of bacteria in parallel to the solid support. Under these circumstances it is foreseeable that the time necessary for the solid support to be covered, at least, with one layer of bacteria decreases as ρ increases. In fact, the roots of ( ) = 0 decrease when ρ increases. It was also predictable that, with the increase of ρ, more overlap layers will form. Therefore, the biofilm will have more different local thicknesses.  The beautiful close-ups of Figs 8 to 11 help to better clarify this example. Fig. 8 shows the functions ( ) with 1 ≤ ≤ 253 ( is invisible).
The maximum capacity of Excel is approximately 254 functions in a single plot. The functions ( ) for = 1, 43, 85, 127, 169, 211 and 253 are equally spaced with respect to n, and they are represented in red for a better visualization of the shape of this whole set of functions. All the others were represented in gray and, because they are so close, they seem to have filled the entire space. However this space represents the close-up only up to just = 200. In Fig. 9, 254 functions are represented to observe the shape of the same set of functions from = 0 to = 1000. They are designated in the legend. The function ( ) is represented in red. This function is very close to the function ( ) in Fig. 8. Therefore, it can be considered that it limits the set of functions already represented in Fig. 8. Those that are also in red, on the right side of the graph, show that for long time an almost periodic growth pattern has been established.
While ( ) slowly decreases (Fig. 7) more and more very small layers are stacked on top of those that already exist, which are also very small. And the areas of solid support already covered grow much more in height than in width. They are like vertical fibbers that increase in length Therefore, in this example, the model can simulate a fibrillar biofilm.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 1 March 2021 doi:10.20944/preprints202012.0673.v2 Fig. 9. Close-up of Fig. 8, near time axis (0 ≤ τ ≤ 1000). Functions : θ0, θ1, θ2, θ3, θ4, θ5, θ6 (θ0 not visible) ... and θ13, θ20, θ27, ... θ13+7n, ... θ1203, θ1210, θ1218, ... θ1210+8n, ... θ1810 and θ1818 (only visible till ≈ θ1014). All the functions are represented in gray, except θ258, θ958, θ965, θ972, θ979, θ986 θ993, θ1000, θ1007 and θ1014, which are in red. Fig. 10 is a close-up taken from Fig. 9 and provides a visualization of the growth pattern before fibrillar growth is reached. Fig. 11 is also a close-up made from Fig. 9 and provides a visualization of the growth pattern of the fibbers, almost only in height. The ( ) functions of this example (Fig. 8 to 11) have an oscillating behaviour. Although positive, the values of the functions with this oscillating behaviour are very small. This detail can be observed due to the performed accuracy and is not surprising in systems with nonlinear dynamics, such this one.  In Fig. 12 (ф = 0 and = 1.5) the entire specific rate growth is carried out in parallel to the solid support and the detachment is strong. It is observed that the system tends, in the cases of = 0.001, 0.01, 0.1 and 0.2, to a well-defined steady state, in which the fraction of area not covered, ( ), is much greater than all the others. In these cases the process of arrival and attachment and the growth parallel to the solid support don´t compensate the strong desorption. Therefore, the biofilm is patchy and consists only of small "hills" or clusters with one layer, two overlapping layers and little else. The close-ups of Fig. 13 show that ( ), ( ) and following functions have also reached steady state. This was already observed in Fig. 12 for = 0.1 ( and ) and = 0.2 ( , and ). In cases with = 1, 5 and 10 the process of arrival and attachment is much more strong.
Consequently ( ) evolves to significantly smaller values.
In this examples also an almost periodic growth pattern is established. The fractions covered with fewer layers are virtually zero and the fractions of greater thickness continue to grow, then reach a maximum and finally decrease because new layers grow over them.
Biomass grows linearly over time due to this almost periodic evolution pattern.
The biofilm covers the entire solid support area, with an external interface on which the statistical distribution of thicknesses will remain unchanged in shape but whose average value increases with time.  Bessel functions of first kind (Case 2) and the borderline solution (Case 3). A recurring algorithm with five terms was built, whose accuracy allows calculations with 12 significant decimal places.
In Case 1, the function ( ) has a root corresponding to the instant when the total coverage of the solid support is reached. It is because full coverage is achieved that it is called continuous biofilm. This happens in all the examples in Fig. 5 and those in Fig. 6 for ≥ 0.1. In Fig. 5, the local thicknesses obtained when the total coverage is reached vary from a single mono layer ( = 0.001) to a thickness distribution between 3 and 7 overlapping layers, with 5 being the thickness corresponding to the largest covered area.
In Fig. 6, for ≥ 0.1, the statistical thickness distributions represent higher values: from an average of approximately 100 overlapping layers ( = 0.1) to 10 ( = 10). Therefore, the structures in Fig. 5 and those in Fig. 6 for ≥ 0.1 vary between the monolayer and a wide variety of local thicknesses that give the outer interface of the biofilm a very irregular shape, with "deep valleys" and "high mountains". However showing a statistical distribution of thicknesses whose shape remains but whose average moves, over time, to higher values. In this case the model is able to translate a traveling wave front into mathematical terms.
We can conclude that, even in its simplest version, the deduced model manages to simulate a wide variety of structures and, additionally, for certain parameter choices, it also goes beyond the initial phase of formation and growth, passing to the phase of total biomass linear growth, over time.

Future work
We can now relate some following deductions that can be done next.
The mathematical description of the linear phase of global biomass growth and accumulation can also be described for continuous biofilm (Case 1).
Once the already described process is started, at = 0, there will come a time when the fraction without any layer of bacteria will cancel out: ( ) = 0. After this instant the solution remains mathematically correct but will not be physically meaningful because will pass, due to the continuity of the analytical solution, to be negative and there are no surfaces whose area has a negative value. In the sequel we designate this solution, valid for 0 ≤ ≤ as "Solution 1".
Consequently, the following instants should no longer be calculated with this solution. But the values of the fractions , with ≥ 1 , obtained at will be the set of initial (positive) values of a new period of time, . Of course, the sum of these new initial values is 1.
In this new time interval a more general analytical solution can be deduced and applied. Now we have several variables whose initial value must be taken into account and not just one. This more general analytical solution guarantees us the analytical continuity in of the model, solving the interval . In the sequel we designate this solution as "Solution 1*" and briefly describe this analytical continuation of the model.
During this second time interval, , it will now be the fraction that will be annulled. When this happens, by the very same reason that there are no surfaces with negative area, we will again consider that the values , ≥ 2, obtained when ( ) = 0, will be the initial values of a third interval, . This third interval, ends when ( ) = 0. And so on: whenever the fraction , with the lowest index (the smallest thickness fraction) cancels out, the following instants should no longer be calculated with the same exact solution since that fraction ceases to exist. Then at that instant, when , with the lowest index cancels, the values (all positive) of all remaining fractions , will be the set of the initial values for a new time period.
Now we can summarize, in greater detail, the sequence of these successive solutions.
-Initial condition for interval 1 (lasting the time and during which will be canceled) That is to say, the successive intervals of order seem more and more to each other. The statistical distribution of the fractions tends to a histogram that, with the proviso that it suffers small variations between the beginning and the end of each interval, remains always the same. This histogram only shifts, over time, to with increasing indexes " " , increasing all of the " " one unit for each " "-interval elapsed. This is how the mathematical translation of a travelling front can be done for continuous biofilm (Case 1). This "Solution 1*" is always the same for all intervals (with ≥ 1) because all the fractions , at the beginning of the interval, have at least one layer of bacteria and this implies that all have growth potential by reproduction , and consequently the global population increases all over the entire surface. This does not happen in the first interval because, during this interval, there is , with a non-zero value ( > 0).
Mathematically this "Solution 1*" is always the same because it starts from a system of differential-difference equations in which the first one explicitly expresses the derivative of following those of , , ... changing only the interval order . This does not influence the deduction that leads to the successive solutions. It is only necessary to properly replace the index as well as the initial conditions that correspond to it.
"Solution 1", applicable only to the first interval, is therefore not exactly just a particular case of "Solution 1*" although this one may be used with any initial condition that we consider convenient. This is so because there are details that distinguish the equation for the derivative of from the equation of the derivative of , for all other intervals, when one put = 1 in this last equation.
The deduction just explained will be next work.
Another future development can be the determination of total biomass evolution over time and also of several variables that characterize the structure.
A mono-layer sufficient to cover the surface unit of the solid support is the biomass unit.
With the evolution over time of all the factions ( ), the evolution over time of the total biomass can be easily And the structural entropy, ... where if one ( ) is zero, the term in the sum is zero.
The greater the structural entropy, the more irregular is the biofilm is and, therefore, the greater the variety of local thicknesses. Formula (A29) is valid for ≥ 1. For = 0, the formula (A19 (0)) must be used for calculate . Considering (A21) we can use, from now on, the coefficients or the . For sake of economy in notation is better to chose, from now on, only the coefficients .
The desired solution for Case 1, in terms of Bessel functions of first kind, ( ), is eq. (80) which can now be found in section 6.2.
For such achievement one only need to apply the early defined formula (81), but inverting, so making explicit

Appendix B. Patchy biofilm solution
Let recall that the solution (117) is only the solution for the named Case 1, corresponding to the continuous biofilm. As signalled before, (117) corresponds to the condition ˂(1 − ).
The "complementary" condition, ˃(1 − ), designated before as Case 2, relates the evolution towards a steady state where the bare fraction don´t never get null. In physical terms this means the formation, at times large enough, of a patchy biofilm.
The solution for Case 2 can be obtained simply from solution (117), being aware that whenever root √1 − − occurs substitution by · √ + − 1 must be done. Consequently, in (117) the variable and the coefficients and are modified as follows.

Appendix C. Frontier between continuous and patchy biofilm solution
The solution for Case 3 can be easily deduced by induction after doing the corresponding simplifications in the initial system of differential-difference master equations, just before any other step. It is so simple that we only indicate the solution, ... which is equal to ( 1( )), as we wanted to demonstrate.