Inverse analysis of material parameters in coupled multi-physics biofilm models

In this article we propose an inverse analysis algorithm to find the best fit of multiple material parameters in different coupled multi-physics biofilm models. We use a nonlinear continuum mechanical approach to model biofilm deformation that occurs in flow cell experiments. The objective function is based on a simple geometrical measurement of the distance of the fluid biofilm interface between model and experiments. A Levenberg-Marquardt algorithm based on finite difference approximation is used as an optimizer. The proposed method uses a moderate to low amount of model evaluations. For a first presentation and evaluation the algorithm is applied and tested on different numerical examples based on generated numerical results and the addition of Gaussian noise. Achieved numerical results show that the proposed method serves well for different physical effects investigated and numerical approaches chosen for the model. Presented examples show the inverse analysis for multiple parameters in biofilm models including fluid-solid interaction effects, poroelasticity, heterogeneous material properties and growth.


Introduction
Microorganisms tend to live in aggregates rather than dispersed single cells and surround themselves with a network of extracellular polymeric substances (EPS) as a survival strategy. This is called the biofilm matrix [1]. Amongst others it provides them with mechanical resistance against external forces acting on them through a surrounding fluid flow. There is a broad variety of biofilm occurrence. Depending on whether they are intentionally grown or unwillingly observed there are opposing interests when it comes to biofilm control in engineered systems. The prerequisite in all cases however is to have good knowledge about their behavior. In engineering good knowledge is represented by the availability of accurately parameterized models that allow the generation of reliable model predictions.
There has been a variety of efforts to estimate relevant material parameters of biofilms. The reader is referred to [2,3] and [4] as they all give an overview of practiced testing methods for parameter estimation for different physical aspects. Therein a variety of mostly intrusive test strategies and results can be found.
As it is still unknown how much influence the environment of the biofilm has on its mechanical properties [5], an effort towards improved in situ testing methods has been made. The long term goal of the methods developed in this work is to develop a nondestructive testing protocol for biofilms to estimate material parameters that can later be used to predict biofilm behavior and develop strategies to control their appearance and growth. Optical coherence tomography (OCT) currently seems to be the prime option to scan the geometrical representation of biofilms in flow cell experiments under variable load [6]. This has led to new insights into biofilm mechanics research [7,8]. The advantage is, that automated growth and test protocols [9] can be developed with this technique and therein the biofilm can be kept in the same environment for the whole test cycle.
First parameter estimation approaches are described in [10]. A common scalar approach for determining stiffness and shear resistance is first presented therein. This type of analysis, if it is applicable, is quick to use and serves with estimates of the values in the relevant order of magnitude. While such approaches in the end might not be sufficient for predictions with suffient accuracy, they can favorably be used to assess good initial guesses for more advanced approaches like the presented inverse analysis algorithm. It has been shown in [8] that fully resolved fluid-solid interaction (FSI) simulations can be used to determine one material parameter from modeling flow cell experiments with linear solid mechanics. This has the advantage that in general no specific shape characteristics of the biofilm are required for the parameter estimation.
The present work proposes an inverse analysis algorithm incorporating different numerical models for non destructive experiments with in situ measurement via OCT. The intuitive approach followed in this work is to model the experimental setup as accurately as possible and then try to find the set of parameters that fits experimental data best, in present case via a Levenberg-Marquardt optimization. In order to best analyze and present the properties of the approach, we create a well defined, i.e. clean environment and hence use artificially generated results via forward numerical analysis and the addition of Gaussian noise as a first step and to particularly isolate important effects arising in this kind of analysis. The honest analysis and discussion of strengths and weaknesses of this type of inverse analysis approach is the main motivation to present this as standalone work.
The focus in this work is specifically set on flow cell experiments, like the ones presented in [7] and [11]. In this type of experiments a biofilm is situated in a transparent channel, the flow cell. The biofilm is grown on the channel floor and supplied with a solution of nutrients in the fluid. Once a substantial patch of biofilm has formed, it can be exposed to varying loads in order to observe its deformation. The mechanical load on the biofilm is controlled exclusively via the fluid volume inflow rate into the flow cell and for growth the composition of solutes in the fluid supply. OCT provides precise measurements of biofilm geometries, but it is limited in the image capturing speed, especially when it comes to three dimensional scans [6]. Empty flow cells were designed as channels with rectangular profile with cross sections in the millimeter range in recent experiments. Exemplarily those had the dimensions of 50 mm × 5 mm × 0.45 mm in [9] or 124 mm × 2 mm × 1 mm in [7].
Experimental results show different shapes of biofilms for different flow rates through flow cells [11]. Even under automated cultivation and measurement, a certain level of variability can be found [9]. The appearing nearly arbitrary shapes of biofilm surfaces pose the question in which way the model results should be compared to experimental results. The variety in biofilm appearance impedes the accessibility to simple or even scalar comparisons for parameter estimation. When it comes to quantifying flow cell experiments, the load on the biofilm is sometimes represented by the fluid shear stress next to the wall of the empty flow cell. It is known [8,12] and quite obvious that this shear stress is not fully representative for the load on the biofilm, as the analyzed biofilm patches in published analyses take up a significant portion of the height of the flow channel. The flow rate itself is a more representative quantity. It can also be represented by average inflow velocity or Reynolds number in the fluid flow through the empty channel. Because of the irregular shapes of biofilms, an interplay of forces and local effects in the fluid flow is affecting biofilm deformation in flow cells significantly. This is addressed in this work by full geometrical resolution of flow cell experiments and a consideration of fluid-biofilm interaction at the interface.
In order to be self contained, this article provides a very brief introduction to physical models applicable to different aspects of biofilm mechanics and a brief summary on numerical approximations. In the next step a method for comparing different surfaces is presented and the optimization used to minimize the deviation of model predictions and experimental results is demonstrated. Then the presented algorithm is tested on different cases and the results are discussed. Finally general remarks and the interpretation of the informative value of results achieved with the presented algorithm are discussed. The article is then brought to a conclusion.

Biofilm modeling
Properties and performance of inverse analysis approaches depend strongly on the model used to represent the experiments that are used to drive the inverse analysis. Results of an inverse analysis can never be better than the physical model incorporated. In the given setting of inverse analysis the numerical model of the experiment is referred to as the forward model. The forward model must include all physical aspects that are relevant. In the same sense all the aspects a model covers should be represented in the experimental data used. The choice of a suitable forward model therefore depends on the type of data used and the type of application for parameters to be determined. The best case is that the model used for modeling the experiment in an inverse analysis is the same as the one, the parameters are aimed to be used with to make valid predictions for biofilm behavior in the future.
In the following we briefly sketch models and methods that we have developed in the past and that are used in this work. Modeling approaches to fluid-solid interaction with scalar transport have been proposed in [13] and have been extended to include biofilm growth in [14]. Porous properties of biofilms and porous flow through them are known for a long time [15]. For solving fluid-poroelasticity interaction (FPI) a novel method was developed in [16] and applied to a kind of finger shaped biofilm example therein. The presented workflow works independent of the material model. We exemplarily choose Saint-Venant Kirchhoff or Neo-Hookean material models in our examples.
The inverse analysis algorithm proposed in this paper is obviously not limited to the kind of models presented here. It can also be applied for other models, like one that describes damage in the sense of detachment [17] or viscoelastic behavior [18].

Physical models
The analyzed type of fluid-biofilm interaction includes significant deformations and potentially rotations of subdomains. Therefore the application of nonlinear kinematics is essential. In addition the observations made include a variety of effects of biofilm physics and the according equations for fluid-poroelasticity interaction, scalar transport of potential nutrients and growth are summarized in the following subsections. For the sake of brevity the presentation is limited to the strong field equations. The respective boundary conditions are implicitly assumed to be well defined. Details can be found in the referenced specific literature.

Fluid field
For the fluid field the general Navier-Stokes equations for incompressible flow of a Newtonian fluid are appropriate. In order to allow for deforming domains, that are essential for fluid-solid interaction, they are given in arbitrary Lagrangean-Eulerian (ALE) formulation as These equations relate fluid velocity v F , ALE convective velocity c F , fluid pressure p F , fluid density ρ F , fluid dynamic viscosity μ F and a body forceb F in the fluid domain F × (0, T ).
Herein the strain rate tensor (v F ) is a short expression for The ALE formulation is one popular approach to model fluid-solid interactions with a moving mesh to ensure the continuity between fluid and solid in case of a moving interface and therein the ALE convective velocity c F is the fluid velocity relative to the moving mesh.

Solid field
For modeling the nonlinear behavior of solids also allowing for large deformations, the general balance of momentum in reference configuration applies. It relates the solid density in reference configuration ρ S 0 , solid displacement d S , deformation gradient F , second Piola-Kirchhoff stress tensor S and the body force in reference configurationb S 0 in the solid domain S 0 × (0, T ). The solid acceleration a S is the second material time derivative of the displacement a S = d 2 d S dt 2 .

Fluid-solid interface
On the fluid-solid interface F,S × (0, T ) balance of tractions h S and the equality of velocities, stemming from a no slip condition between fluid and solid, need to hold. The solid velocity v S is the first material time derivative of the displacement v S = dd S dt .

Scalar transport
A scalar transport model in the coupled fluid-solid model is used for nutrient distribution during the flow cell experiment. The solution of this type of systems has been presented in [13] and [14]. The scalar transport equation is stated for the fluid (5a) and the solid (5b) domain.
For the domains the concentration of the solute is described as K in fluid and solid domains. The respective diffusion coefficients are described as D K , with K ∈ F, S being the index for quantities in one of the phases. R S is the reaction rate. In this work a Monod kinetic relates to nutrient consumption of the biofilm which therefore only appears in the solid domain.
This kinetic includes two coefficients, which are the reaction rate K R 1 and the half saturation K R 2 . The negative normal flux on the interface is computed as for phases K with the respective unit normal n K . On the fluid-solid interface the concentrations and fluxes of the phases K must be equal.
In this formulation (8b) the fluxes on the interface must be related to the same normal direction on the interface. The scalar transport is only one way coupled to the fluidsolid interaction as the deformation of the solid and the fluid velocity will influence the concentration solution, but the concentration does not really influence solid or fluid field solutions.

Surface growth
We introduce growth or erosion according to [14]. It is assumed that growth or erosion is localized on the biofilm surface and influenced by nutrient flux and surface tractions.
With n S being the outward pointing normal on the biofilm surface and t S i two respective orthogonal unit tangents. This is a relatively simple phenomenological model wherein stress inhibits growth or erodes the biofilm [14] and the nutrient flux into the biofilm domain is the cause for growth. The normal flux contributes to growth with the factor K g 1 and the erosion is determined with the factors K g 2 due to the normal stress component and K g 3 due to the tangential stress components. t g describes the timespan the biofilm is exposed to given growth conditions and d S g the resulting displacements of the surface because of growth.

Poroelasticity field
The poroelasticity model, that is used, relies on the assumption of a homogenized mixture between a fluid phase with Darcy flow and a solid structure. Porosity φ acts as volume ratio of void that is filled with a fluid phase. This is well described in [19] and a fully coupled numerical model for the field equations was developed in [20]. The coupled system of equations for the poroelastic mixture is derived as In (10) the indices (•) P S for the porous structure phase, also called skeleton and (•) P F for the fluid phase are used. J = det F is generally known as the Jacobian determinant being the determinant of the deformation gradient F . The time derivates ∂(•) ∂t X P in (10) are evaluated in the reference configuration of the porous domain with the reference coordinate X P held constant. The skeleton velocity v P S and skeleton acceleration a P S are the first and second material time derivatives of the skeleton displacement d P S . With the presented formulation the porosity needs an initial value φ 0 for each material point and then changes according to the skeleton displacement field of the fully coupled model during deformation.ρ 0 P S = (1 − φ 0 )ρ P S 0 represents the macroscopic averaged initial density of the skeleton and k = (J ) −1 F · K · (F ) T is the spatial permeability computed from the permeability K in reference configuration, which is determined with the Kozeny-Carman formula.
The porous composite strain energy function is Therein ψ P,skel is the contribution from the skeleton, ψ P,vol accounts for the volume change due to changing fluid pressure and the penalty part ψ P,pen guarantees positive porosity. Those are as in [20] with the initial porosity φ 0 . This formulation allows usage of any hyperelastic material model for the skeleton part of the strain energy function ψ P,skel . We choose the scaling parameters as penalty parameter η = 1.0 and bulk modulus κ = 100 accoring to [20] for the presented example.

Fluid-poroelastic solid interaction
A consistent approach to tackle the interaction of a Newtonian fluid and a poroelastic solid is presented in [16]. The method presented therein is directly applied in this work. On the interface must hold. The conditions describe a balance of tractions between the porelastic mixture and the pure fluid (14a), equality of fluid pressure in the poroelastic and fluid domain (14b), the continuity equation for the normal fluid flow (14c) and the so called Beavers-Joseph conditions [21] for the coupling of the tangential components i = 1, 2 in directions t F i of the fluid velocities (14d). Therein the interface permeability κ is used. The Beavers-Joseph constants α BJ , β BJ regulate this tangential velocity dependency in (14d). They are both chosen to α BJ = β BJ = 1 in the presented example.

Numerical approximation
For all numerical models a nonlinear finite element method (FEM) based approach is used. For time integration a one-step-theta approach is used. For the pure FSI examples we use a monolithic arbitrary Langrangean-Eulerian (ALE) approach just as in [13] and [14]. For the mesh deformation the ALE fields can be treated as a quasi-elastostatic pseudo solid. Monolithic methods are preferable for FSI problems in many biological applications as they might contain fields with similar density with soft solids, represented here by low Young's modulus [22]. ALE Methods can be problematic when it comes to a change in the topology or large mesh displacements. An alternative approach to overcome difficulties associated with these cases are fixed grid methods. For the fluid-poroelasticity interaction examples presented, a CutFEM based approach [16] is used. For FSI problems an approach based on the so called CutFEM has been developed in recent years [23]. It is capable to solve FSI problems with a fixed fluid grid. For a fixed fluid grid the fluid equation (1a) must be replaced by the Euler formulation with zero grid velocity. Thus, the ALE convective velocity c F is replaced by the fluid velocity v F for the fixed grid. The idea is to cut out the parts of the fluid mesh that are covered by the solid and solve the fluid field on the remaining discretization including partially cut elements at the interface. To sustain a proper fluid mesh with uncut elements in the interface vicinity also a hybrid method of ALE and CutFEM was developed in [24]. CutFEM based FSI approaches enable the treatment of cases when it comes to a change in topology like for partial or full detachment or on the other side self contact [25,26].
For the growth algorithm there is the need for multiple time scales as the FSI dynamics acts in the range of seconds and the growth processes take place in the range of days. For this multi-scale approach in time a quasi steady or periodic state for the fluid-solid-scalar interaction is reached with smaller time steps. Based on the FSI and scalar transport solution the nutrient flux and interface tractions are evaluated and surface growth is applied with the much longer growth time step. This procedure is then advanced until the full growth time period is reached. After the growth step an ALE relaxation of the mesh is necessary for the domains on both sides of the interface fluid and biofilm, to distribute the displacements smoothly on the whole domains and thereby reduce mesh distortion [14].

Surface distance measure
We set inverse analysis as a special optimization problem and in the application of optimizers the quantity that is minimized is called objective function [27]. In inverse analysis the objective function somehow describes the difference between experimental observation and forward model output. In this context it turns out that the measure that is used to quantify this difference is a key question in inverse analysis. The result of an inverse analysis always depends on the approach used for this measure. Hence, detailed information about the measurement approach need to be combined with the achieved results for presentation and interpretation. This also makes it obvious that the selection or design of a suitable measure is crucial and must be well considered. One important contribution of this work is to propose a simple geometric measurement for the objective function in this kind of experimental settings for biofilm parameter estimation problems, that is suitable for any optimizer.
As OCT measurements of experiments only contain information if a point in space is likely covered by biofilm or not, there is no pointwise displacement information available. Comparable pointwise displacements from the observed experiment would need to be somehow computed first. But in order to do so, some additional assumptions would need to be introduced, which in turn spoil or bias the outcome. In addition those assumptions -being more or less physical -might even substantially complicate the inverse problem or they might point to "unphysical" scenarios. Because of this we argue that it is not the best idea to compare the forward model evaluations to a field of selected point displacements, which are themselves the result of a postprocessing operation, but rather refer to some primary information, which, in the case of flow cell experiments and OCT, is the surface shape. Evaluating this information is performed as depicted in Fig. 1.
Given an observed result of an OCT measurement the first step is to determine a representation of the fluid-biofilm interface obs , by some sort of image segmentation. Image segmentation is a topic on its own and not the focus of this work. For the purpose of this paper, it is enough to assume that the data is already suitably segmented, which can also mean a segmentation done by hand. On the observed interface obs the analyst is to choose significant points, meaning points where the interface underwent significant displacements during the experiment. These measurement points are depicted as crosses in Fig. 1. The distribution and number of measurement points is up to the choice of the analyst as it depends on many aspects. The number should at least be greater than the number of parameters, that are to be determined in the inverse analysis for the optimizer Fig. 1 Prototype sketch for pointwise distances d i mp from measurement points on the deformed interface obs , observed from experiment, to M resulting from forward model evaluation, in given measurement directions depicted as rays to work well. The actual choice of measurement points should in our impression be made towards regions, where the validity of all model assumptions is trusted the most. That means it should favor points away from potentially uncertain boundary conditions and potentially uncertain flow conditions in the channel and towards regions, where the spatial resolution of OCT and the quality of the captured images is trusted the best. Secondly, given the forward model and the parameters analyzed, the measurement points should be chosen in a way that all the different parameters can show significant effect in the resulting distances.
For every point selected, an individual search direction for the intersection with the result for the displaced fluid-biofilm interface from the forward model evaluations M must be decided. These are depicted as rays in Fig. 1. A general recommendation is the normal direction. Nevertheless again the quality of the image decides if a confident guess for that normal can be made. Given the fact, that OCT scans are generated from above the experiment, the vertical direction is another reasonable choice.
With both the measurement points and the associated directions at hand the thereby defined rays must be intersected with the deformed interface resulting from forward model evaluations M . As the resulting distance for every measurement point d i mp is based solely on a geometrical measure, it has no inherently conclusive sign. Therefore it is defined positive if the intersection point is on the side of the biofilm with respect to obs and negative if it lies towards the outside. In the case of multiple intersection points the lowest resulting distance is used.
The presented choice of comparative measure for the interfaces including both measurement point location and measurement direction circumvents known drawbacks of the closest point projection described in [28]. With the proposed method the distance measurements are uniquely defined as the search direction is predefined. This can be very useful as, if two or more candidates for the closest point on M to a measurement point exist, a gradient based optimization with a finite difference approximation as the one presented, can be heavily deteriorated. The capturing of irrelevant shape characteristics must be prevented by a good choice of measurements points by the user.
Overall the presented measurement method is considered rather hands-on, as the observed interface location must only be determined for the measurement points. For our type of problems, this is also a clear advantage compared to the usage of global surface comparisons as for example the ones presented in [29] and [30]. For global approaches for surface comparison a full representation of the surface must be available and therefore must be constructed from the data. A global measurement approach also poses higher demands on the image segmentation than the presented method. Nevertheless in the case of optimal data acquisition and the assumption that the experiment is modeled optimally, the presented measurement is also fully automatable using factual normals e.g. for every triangle of a triangulation of the observed deformed biofilm surface obs .

Levenberg-Marquardt optimization
For the minimization of the objective function defined by a suitable comparison of observed experiments and a forward model we use a Levenberg-Marqaurdt approach for optimization. Levenberg-Marquardt optimizers go back to the works of Levenberg [31] and Marquardt [32]. A good overview of the actual algorithm is shown in [33]. A Levenberg-Marquardt optimizer is in general a deterministic method. As a gradient based method it represents a local optimizer and is applicable for inverse analysis if the dimension of the inverse problem is low enough and the initial guess is good enough. In the selected numerical examples we will shed some light on the applicability for different numbers of parameters and initial guesses for our target applications. In the past we have already succesfully applied such algorithms for identification of constitutive laws and parameters of hyper-and visco-elastic biomechanical problems (see e.g. for problems with single type of experiments [34,35] and also the combination of different experiments on the same specimen [36]).
In order to have a rather self contained paper we will briefly sketch the algorithm in the following. The Levenberg-Marquardt method is used to minimize a least squares objective function with forward model results M (x) j at potentially different times for n x unknown parameters in the parameter vector x and n r observed experimental measurements (y obs ) j . The algorithm uses the regularization parameter μ and is started from an initial guess x 0 , μ 0 . The core algorithm is to iterate the update rule for the parameter vector with the parameter step computed by until predefined convergence criteria are met. The iteration index k is omitted in J k , μ k and r k from here, as it is clear that terms should be computed exclusively at current step k. We propose that the vector of punctual distances d i mp , potentially also collected over different time steps, between the forward model result surface and the observed experimental surface measured in the way presented should be directly used as residuals and arranged in the residual vector r of length n r . In the least squares formulation of the objective function (16) The Levenberg-Marquardt method makes good use of the vector shape of the residual for finding the parameters for the next step. The algorithm uses the partial derivatives of the residual components with respect to the parameters as the Jacobian In absence of actual gradient information, the Jacobian is approximated by finite differences. For that, n x + 1 simulations per iteration are necessary to be computed. One model evaluation is conducted with the current parameter set x (= x k ) . n x further model evaluations are computed with the i th parameter perturbed tõ Resulting n x perturbations of the parameter are written as vectors x i . The approximations of the partial derivatives are then arranged into the approximation of the Jacobian To assess how much the current step has improved the result, the gradient based error err k grad J T r 2 k := err k grad (24) can be used. In order to simply evaluate how close the current model solution is to the experiment, the residual error err k res can be computed as The regularization parameter is updated according to only if the current step is closer than the previous step.
Treating the regularization parameter in this adjusting way is the only form of adaptivity in the algorithm. In general it helps to find a suitable step size by reducing the regularization especially close to the optimum.
So far the presented Levenberg-Marquardt method shares the property with its family of optimizers to be unbounded. No information about the validity of parameters is introduced so far. Often material parameters come with a valid interval, or constraints, that need to be fulfilled. The pure algorithm presented so far is not constrained, so potentially if the Jacobian indicates further decrease of the residual into a given direction, the algorithm suggests parameters x k+1 , that cannot be used in the model. On top of the Levenberg-Marquardt algorithm we want to be able to set constraints on every parameter.
To achieve this property an additional check is introduced. If any suggested parameter in x k+1 is out of bounds, the step is declined, the regularization parameter is doubled μ k = 2μ k and a new step x k+1 is proposed. To avoid useless model evaluations the algorithm is terminated if μ k grows unreasonably high μ k > μ 0 · 10 6 . In that case no parameter result can be found. The algorithm is terminated if a certain convergence criterion is met or if a maximum number of iterations is n max reached For real experimental data, comparison to grad is the better suited criterion, because the residual error in the data is unknown. This way iterations are stopped, when the gradient does not indicate significant decent in the residual measure. As err grad has by definition (24) no unique physical unit as it depends on different parameters in x, physical units are omitted for this quantity. Overall the algorithm is deterministic and a local optimizer. So, if more than one local optimum exists within or outside of the bounds, there is no guarantee to find a global optimum even for a continuous and bounded problem.

Numerical examples
In the following numerical examples we want to show that given inverse problems in biofilm physics are well solvable with the presented measure for the similarity of surfaces and the Levenberg-Marquardt approach. The general setup is representing a prototype flow cell experiment, wherein a solid representing the biofilm is exposed to a certain volume flow rate from the left and is therefore deformed towards the right. As already stated above, the performance of an inverse analysis depends on a number of things beside the method itself, like the question how well the numerical model represents reality in the experimental setup. Hence in order to get a better impression of the quality of a specific method for a certain type of application, it is advantageous to test such methods on "clean data" first. A common approach to generate such data is to use the numerical model in a forward analysis with some chosen parameters and potentially add some noise to the results, in order to generate artificial experimental or measurement data to be used in the following inverse analyses.
The numerical examples were computed using the referenced methods implemented in the inhouse multi-physics C++ code BACI [37] and a tailored python framework QUEENS [38]. The presented inverse analysis algorithm was newly implemented in QUEENS during this work. QUEENS is used to run and manage forward model evaluation and conduct the inverse analysis. The intersections of the rays representing the measurement directions and the mesh based forward model results have been found with the python package for vtk [39].
Although The flow cell experiments are modeled to last several seconds. The focus is set on the quasi static case, meaning that the biofilm is free of oscillatory or inertia effects in the observed reference results from forward model evaluation. For that, the inflow rate is applied smoothly with a cosine based function in multiple steps of an increasing volume inflow and then held constant. The inflow is assumed to be the result of laminar channel flow and therefore chosen to have a parabolic profile. This is how the inflow boundary condition can be reduced to a single quantity, the volume inflow rate. On the right hand side boundary a horizontal outflow is enforced, to reflect, that it is no free outflow, but the channel continues downstream from the modeled region.

Homogeneous biofilm
As a first and most simple example the presented inverse analysis algorithm is performed for a fully homogeneous solid model of a biofilm interacting with the fluid flow. The shape of the biofilm domain in the simulated model is arbitrary and inspired by experiments shown in [7,8]. The modeled biofilm patch is held in place on the channel floor with a no-slip condition on its lower, straight boundary. A parabolic fluid flow profile with flow rate 100 mm 2 /s from the left boundary of the purely two-dimensional channel with 2 mm × 1 mm is ramped up for 10 s. After 15 s a quasi steady deformation state was observed and thus the deformed state is regarded at that time.
For the presented FSI examples a Saint-Venant-Kirchhoff material is used like in other biofilm related works [17,40]. Its behavior is governed by two parameters, namely Young's modulus E and the Poisson's ratio ν. This material is linear in the Green-Langrange strains and second Piola-Kirchhoff stresses and for cases with small deformations can be related to the classical Hooke constitutive law, which is standard in linear continuum mechanics and also used in other biofilm mechanics studies [8]. The biofilm is obviously threedimensional and it is assumed that its behavior is not changing much in the out of plane direction, which is reflected by using a plain strain assumption for the reduction to two dimensions.
First, the reference result with parameters E = 400 Pa and ν = 0.3 in the biofilm model is computed and the field solutions shown in Fig. 2 are obtained. The biofilm domain fills the volume between the channel floor and its upper surface, that is the fluid-biofilm interface. In Fig. 2b additionally the interface tractions acting on the biofilm, resulting from the fluid-solid coupling, are displayed as arrows. The simplistic assumption commonly used in biofilm mechanics of a constant tangential force on the whole fluid biofilm interface would be distinctively inaccurate in this example, as it can be observed that the interface    With the given prerequisites the inverse analysis of the two material parameters Young's modulus E and the Poisson's ratio ν is conducted with different initial guesses listed in Table 1. The cases with negative Poisson's ratio are included as there is some speculation in the literature wether biofilms are so called auxetic materials. Resulting search paths in the parameter space are shown in Fig. 4 with the associated colors, which are consequently used throughout this example. In Fig. 4 and all following plots each marker represents one Levenberg-Marquardt iteration.
To assess the impact of noise in the data on the inverse analysis result, different scales of normally distributed (Gaussian) noise with standard deviations σ of 10 −4 mm, 10 −3 mm and 10 −2 mm have been added to the measured points representing the experimental data. As compared to the displacement field of the reference result shown in Fig. 2a with a maximum displacement Magnitude ≈ 5.3 · 10 −2 mm, that is about a fifth of that or lower. For all data sets the algorithm has been run for the same initial guesses. The results are summarized in Table 2 for statistics over all algorithm runs with all different initial guesses for the individual noise levels and the noise free case. With increasing noise on the data the remaining residual error err res increases. Means of the parameter results drift away from the ones used for the reference forward model evaluation and the respective standard deviations in the parameter results increase.
The actual resolution of OCT measurements is in the range of μm [6,9]. In this simplest of the presented examples and this primitive uncertainty estimation in Table 2 it appears that the expected accuracy for the inverse analysis with this level of noise reaches at best the first two digits of the parameter results.
From the residual and gradient based errors over the iterations, plotted in Fig. 5, it can be seen that those quantities decrease steeply for the noise free case, when the parameters come close to the local optimum. Only the step size used in (21) for the finite difference  (c) Fig. 6 Development of a residual error err res , b gradient based error err grad and c regularization parameter μ for inverse analysis runs with different initial guesses over iterations, with Gaussian noise with standard deviation σ = 10 −4 mm added to the measurements approximation limits the accuracy in this setting. For the noisy data it can be seen in Fig.  6 that there is a clear limit to the achievable residual errors err res and some diffuse barrier for the gradient based errors err grad already for the slightest noise. To show this effect the convergence criteria were intentionally not adapted, although it is obvious from Fig.  7 for σ = 10 −3 mm that grad = 10 −6 would have been a good choice, because there is a distinct level for err grad that cannot be reached even with many more iterations due to a low convergence criterion grad . For the data set with largest noise level shown in Fig. 8 the gradient based error couldn't be reduced to a tight convergence criterion of grad = 10 −8 , but an adapted value of grad = 10 −6 did lead to convergence for all initial guesses. It can be further observed that the numbers of iterations until a feasible level of err grad is reached does not vary much for the same initial guesses. From the third column of graphs in Figs. 5, 6, 7, 8, that show the regularization parameter, it can be observed, that as expected the regularization adapts towards low values close to the result, to accelerate the progress of the iterations for all inverse analysis runs. The path in parameter space and therefore the assumed overall shape of the residual error err res in parameter space does not change significantly for higher noise levels, as seen in Fig. 4, although the level of remaining residual error changes in orders of magnitude. Nevertheless, the reached optimum changes significantly for the Young's modulus E = 506.2 Pa instead of 400 Pa for the maximal noise level, but is rather insensitive for the Poisson's ratio. It is obvious and also shown in Fig. 8 that for the data set with the highest noise the level of remaining error is also the highest, but the algorithm still converges to the same local optimum for all chosen initial guesses. Mind that error plots are all presented in logarithmic scaling. So although there is far less progress in the residual error for noisy data, the algorithm finds the shifted optimum repeatedly for all initial guesses.
Looking at the "olive" green and gray lines in Figs. 4 and 5 it can also be observed, that the residual error is low and rather immobile in a local vicinity of the starting points with low Young's moduli and significantly negative Poisson's ratios. Especially for the initial guess of E = 400 Pa, ν = −0.9, i.e. the "olive" green line, convergence is inhibited by this indifferent shape of the residual error for low Young's modulus E and negative Poisson's ratios ν < 0. But there is a physical explanation for this as with this type of interface location based measurement, quite similar surfaces can be obtained via large bending due to a low value for E or by lateral deformation resulting from negative ν.

Heterogeneous biofilm
As indicated in the literature [12,41,42] biofilm material properties depend on growth regimes, age and on induced flow rates. This implicates a layer like structure that a biofilm might develop under varying conditions. An arbitrary showcase example model, with subdomains as depicted in Fig. 9, is used to show that the presented method is applicable to determine different material parameters for different subdomains.
The material parameters in this reference solution are: E 1 = 500 Pa, ν 1 = 0.2, E 2 = 200 Pa, ν 2 = 0.1, E 3 = 1000 Pa, ν 3 = 0.3 in a Saint-Venant-Kirchhoff material model. The geometry used for the homogeneous biofilm model is reused here. The channel geometry and fluid volume inflow are controlled in the same way. The deformation of the biofilm under given load due to the interacting fluid forces can be seen in Fig. 10.
Measurement points and directions are also introduced in Fig. 10 on the deformed geometry. It must be taken care to measure all the influences of all subdomains and therefore two measurement points were chosen on the interface of the stiffer footing layer. Measurement was conducted normal to the observed interface.
At first it is verified that with the noise free artificial measurement data, the method allows recovering the correct material parameters. For that a short summary of algorithm runs with different initial guesses listed in Table 3 is plotted in Fig. 11 with the respective colors. Resulting search paths for number of parameters n x > 2 can no longer be interpreted visually with respect to the response surface in err res . Therefore search paths are plotted individually for the parameters. For these plots one color codes one inverse analysis run. In Fig. 11 it can be observed that the inverse analysis is able to find the reference parameters for noise free data for different initial guesses. It was observed that this problem converges faster if Young's modulus E is underestimated in the initial guess. These algorithm runs are the ones with the highest number of parameters presented. It is obvious in the plots, that the search paths for all parameters are interdependent. This is  also expected from the algorithm as for every iteration a finite difference approximation of the partial derivatives in all parameters is used to find the next step. The great increase in convergence speed seen in Fig. 11a below err res = 10 −4 mm is a hint, that a very local optimum is found, as also the parameters do not change much for those respective last iterations.
Real OCT resolution is in the range of μm [6,9], so the following examples will be run after Gaussian noise with standard deviation 10 −3 mm was added to the generated artificial measurement data.
Remark 1 (Estimation of initial guess) For this setting several arbitrary initial guesses did not lead to convergence of the method for the noisy data. That is why the problem is first run with a reduced set of parameters. To do so and set an application oriented scenario, where the heterogeneous character is unknown and the individual parameters are unknown, the domain is assumed to be homogeneous and the material parameters that fit that assumption are searched for. This also helps to loosen the strong one-toone relationship between generated data and forward model evaluations. The results are displayed in Fig. 12.
In Fig. 13 it shows that the convergence criterion of err grad < 10 −9 mm could not be met. So the result must somehow be concluded from the iterations. It is recommended to judge by an adjusted, but still objective convergence criterion. Therefore the data criterion is adjusted to err grad < 10 −6 mm. With adjusted convergence criterion the averaged result for two different initial pairs of values, listed in Table 4, is E = 323 Pa and ν = 0.0656. For E that is well in between the parameters used for the reference data and for ν that is below all the original values. From this simple numerical experiment we can already conclude that the achievable level of gradient based error err grad cannot be predicted and hence the convergence criteria must be tuned to the data used. Nevertheless the result achieved with the algorithm is conclusive even if the algorithm did not converge under too strict criteria.  In the next step it can be assumed that the domain is in fact heterogeneous and ν is equal for all subdomains. To set this example ν 1 = ν 2 = ν 3 = 0.1 is rounded from the result with the assumption of a homogeneous domain in Remark 1 and E = 300 Pa is picked as an initial guess for an inverse analysis for E in the three subdomains. This results in a distribution of E 1 = 497 Pa, E 2 = 226 Pa, E 3 = 510 Pa. The correct tendency in stiffness in the subdomains is apparent. Only in the footing layer the result is far off the reference value. Most likely the influence of the stiffness of the footing layer is not conclusive enough towards the shape based surface comparison. The remaining residual error is err res = 8.98 · 10 −4 mm. That is lower than the solution with the homogeneous approach and even lower than the residual for a model evaluation with the parameters used for the noise free reference result. This means the added noise has altered the data in a way, that it does no longer represent the reference result in the position of the optimum. If as a further step this result is used as an initial guess for a new optimization for all six Parameters, the residual error can be lowered to err res = 7.43 · 10 −4 mm with the result E 1 = 273 Pa, ν 1 = − 0.60, E 2 = 192 Pa, ν 2 = − 0.10, E 3 = 491 Pa, ν 3 = 0.39. It appears this type of problem and the measurement only via interface deformation favors the assumption, that the material is auxetic, i.e. ν < 0.0. It appears that in this numeric experiment the field of residual error has flipped and the optimum has shifted towards softer, auxetic materials. It is a further conclusion that optimizing for E only is more robust, than the combination of both E and ν at once because a negative ν can compensate a E that is too low in the surface measure. That means that lateral expansion will fill the gap to the optimum for too much bending. Auxetic materials are very rare and mostly occur in specially designed materials with very unique microstructures. As long as that is not proven for the material of interest it is a valid assumption that ν is positive. If nevertheless an optimization result is obtained, that does not seem trustworthy, e.g. negative Poisson's ratio, different initial guesses or even different type of optimizers should be applied and resulting optima compared by their respective residual value err res and plausibility. If the noise in the measurement was too high, or the data not conclusive enough, an increase in data for the specimen might be necessary.
It is known that inverse analysis not only depends on the physical problem at hand as well as on the forward models used, but also on the type and amount of measurement information. For the above type of problems, simply more and other measurement input would be needed in order to allow for identification of all values for Young's modulus and Poisson ratio at the same time. More data could be gathered in form of different snapshots in the same experiment with different load levels or an increase in measurement points. Basically that can be measurements for different settings of the experiments with different load (inflow volume rate) states and developments or different measurement points based on the presented shape comparison. If available, that can also be different measurements from potentially different imaging techniques.

Two-phase poroelastic biofilm
The porous nature of biofilms is well documented. Measurements from OCT scans allow an estimation of biofilm porosity [6]. In the following the attempt to determine porosity via inverse analysis will be presented. The reference result that serves as dummy experiment is set up in a similar manner to previous FSI examples. The biofilm and channel geometry are the same. We switch to the quasi two-dimensional setting with thickness 0.01 mm and use the same inflow rate 100 mm 2 /s · 0.01 mm = 1 mm 3 /s over the height of the left boundary. All displacements and velocities in thickness direction are restricted to zero. Further parameters to the fluid-poroelasticity interaction are arbitrarily chosen as permeability K = 10 −4 mm 2 and Poisson's ratio ν = 0.3. The fluid phase in the poroelastic domain also has the properties of water. For the reference result the parameters E = 300 Pa and a homogeneous initial porosity of φ 0 = 0.25 are used in a Neo-Hookean material model. As we are using a fully coupled two-phase poroelastic model, porosity changes   15 Results for a fluid velocity and b fluid pressure of reference simulation illustrated on the deformed geometry both for "external" flow field and fluid inside the porous medium due to deformations caused by interaction forces from the external flow field but also due to pressure within the porous medium itself. The solution for the velocity magnitude, displacement magnitude and porosity of the biofilm are shown in Fig. 14. It is observed that the biofilm bends with the flow to the right. This is depicted in Fig. 16 where the edges of the initial and deformed geometry ale plotted on top of each other. The porosity opens up in the stretched upstream part and is reduced in compressed downstream regions (see Fig. 14b). The results for velocity and pressure solution of the fluid, inside and outside of the biofilm, are depicted in Fig. 15.
The measurement points used for the inverse analysis are displayed in Fig. 16 on the deformed geometry. They are chosen where the most significant deformation of the interface shape is expected. Knowing, that the varying porosity is coupled to the biofilm deformation one measurement point is chosen on the lower right bump of the geometry. On that basis the inverse analysis is conducted for different initial guesses listed in Table 5. The paths in parameter space are shown in Fig. 17, wherein each marker represents one Levenberg-Marquardt iteration. For two initial guesses, namely the blue and green line, the algorithm converges to the reference parameters and for one initial guess, displayed in orange, the algorithm had to be terminated without result.  It appears that the measured deformation of the interface is not fully conclusive towards the biofilm material porosity, as higher porosity, due to the interplay between porosity and Young's modulus, also lowers the effective stiffness of the porous medium. And as the porosity is naturally bounded φ ∈ [0, 1], the algorithm tends towards the upper bound. Since the gradient indicates further improvement towards this bound, the algorithm gets stuck in the applied upper bound of φ max = 0.9 and the upper limit for the adaptive regularization parameter terminates iterations. Nevertheless, if the initial guess for Young's modulus is good enough, the optimum can still be found, although it takes many steps. It can be observed in Fig. 18, that also the convergence speed depends strongly on the path in parameter space and therefore obviously also on the quality of the initial guess. Along the first (blue) path the analysis converges in nine steps, whereas the one with a higher starting value for the porosity (green) takes 27.
It would also be desirable to include more parameters into the inverse analysis and, for example, to additionally optimize for permeability K or Poisson's ratio ν in the same algorithm. However this short example is only meant to serve as a proof of concept for the suggested approach. In case many parameters and their interplay need to be considered, it definitely makes sense to also include sensitivity analysis and probably also consider probabilistic based (inverse) analysis approaches.

Surface growth of biofilm
This last example demonstrates that the method is also applicable to identify parameters in growth models as the one developed and used in [14]. The inflow rate 0.1 mm 2 /s · 0.01 mm = 10 −3 mm 3 /s, that is induced from the left side, is much lower in this example as it is applied over a long time period providing growth conditions for the biofilm. Growth processes take place on a different time scale than dynamic FSI and this is accounted for in the temporal multi-scale approach detailed in [14]. The geometry of the problem is chosen as the one presented in [14] to sustain comparability. The flow channel has the dimensions 0.6 mm × 0.3 mm × 0.01 mm in a quasi twodimensional model with no displacement or velocity in thickness direction. The biofilm is represented by a finger like structure of 0.04 mm width and 0.1 mm height that ends up with a semi circular tip. The inflow velocity profile is parabolic according to the other examples. The FSI time period is 5 s and the fluid inflow rate is increased smoothly. The growth time period is one day. The material parameters for the Saint-Venant-Kirchhoff material are E = 100 Pa and ν = 0.3. They are assumed to be known from a deformation experiment. Like in [14] the concentration of the scalar species at the inflow is chosen in the range of oxygen dissolved in water as F in = 2.5 · 10 −11 mol/mm 3 . The reaction rates in (6) are assumed as K R 1 = 3.0 · 10 −11 mol/(mm 3 s), K R 2 = 3.0 · 10 −12 mol/mm 3 . The Diffusion coefficient for both phases is D S = D F = 2.5 · 10 −3 mm 2 /s. The scalar in the scalar transport problem represents a dummy nutrient for the biofilm and will be referred to as such in the following.
Model equation (9) shows that the used growth model depends on three different parameters, namely K g 1 as the factor for growth on the domain boundaries scaling with the actual nutrient flux, K g 2 as the factor for inhibition of growth due to normal stresses and K g 3 as the factor for inhibition of growth due to shear stresses. They all appear linearly in the chosen surface growth model. For the reference result the parameters K g 1 = 6 · 10 4 mm 3 /mol, K g 2 = 5 · 10 −2 mm 2 s/g and K g 3 = 8 · 10 −2 mm 2 s/g were used. The reference solution is   Fig. 19 for the fluid velocity and pressure and displacements due to hyperelastic deformation and growth and for the scalar transport species concentration in Fig. 20. The edges of the deformed (bended and grown) biofilm and the initial shape are plotted in Fig.  21.
The changes in biofilm geometry due to displacements and due to growth range in the same order of magnitude. In Fig. 19b surface growth is displayed at the interface along with the solid ALE field on the biofilm domain. This gives a more intuitive overview of the growth deformation, although the ALE field has no physical meaning and the plotted growth stems from a pure surface growth model. The solid ALE displacement field is a reaction to the displacements on the interface induced by growth (9). The displacements are prescribed to the interface nodes of the solid ALE field which is subsequently solved. Figure 19b shows the displacement solution of that solid ALE field which includes the growth displacements on the interface. The nutrient that is transported through the flow is consumed in the biofilm domain according to the reaction rate in (6). This produces a gradient over the interface and therefore nutrient flux, that leads to growth. On the upstream side the growth because of nutrient supply and the inhibition of growth because of higher tractions are more balanced, whereas on the downstream side, where the fluid   induces lower tractions, a larger growth is observed at the interface, although the nutrient flux is lower. Over the finger tip the erosive effects of the traction can be observed. The interface deformations are measured in the points displayed in Fig. 21 on the fully deformed geometry. The distribution of the measurement points is based on the anticipated different regimes -for growth and FSI -on the upstream, downstream and tip side. There must be points in regions with large growth resulting from high nutrient flux or low interface tractions, in regions with high tangential components of the interface tractions and regions with high normal components of the interface tractions.
Initial guesses used for the inverse analysis of the growth parameters are listed in Table  6 and results obtained are plotted in Fig. 22. The errors plotted in Fig. 23 decline over the optimization iterations. It is observed that the search path depends on all growth parameters at once. That means that the residual error err res , as the norm over surface distances, does not measure the three contributions of the parameters to growth and erosion independently in this analysis. This problem setting allows to have initial guesses further away from the optimum, as long as growth is significantly smaller than the dimensions of the biofilm domain. The parameters from the reference data are the result of the inverse analysis for all analyzed initial guesses.

Discussion
Beside the discussion on presented examples also some general remarks and discussions regarding general aspects of the approach that cannot be shown in examples are added for the completeness of the presentation. Further, anticipated general aspects regarding the application of the method are summarized.

Significance of parameters
We have looked at coupled models, that displayed different shapes of residual errors over model iterations in parameter space with the given measure of the biofilm surface shape.
In several combinations, compensation effects in the parameters occurred (e.g. E and ν or E and φ). The presented approach can only be used if the information to all parameters analyzed are actually at least somehow represented in the data and also show effect in the biofilm surface shape in at least partially independent patterns. For presented physical biofilm models the key parameters could nevertheless be determined in example inverse analyses.

Model response
Presented method is not designed to explore the full parameter space, which would be interesting for the global view on the plausibility of combinations of parameters in the whole parameters space for given data. For exploring the residual error on an interval in the parameter space there is a great variety of random or quasi random sampling methods.
To efficiently get an global estimation of the objective function in the full parameter space regression methods like for example Gaussian processes [43] on respective sample results can be applied. The advantage of the presented method is, that it is not necessary to explore the whole parameter space but to find a local minimum with a limited number of iterations and therefore limited amount of forward model evaluations.

Deterministic character
It should be emphasized that the Levenberg-Marquardt optimization is a deterministic approach and there is a risk that the numerical model used for an analysis is not stable along the full search path and especially in the vicinity of the optimum. If the forward model cannot be evaluated and therefore no measurable results retrieved, the finite differences cannot be computed and the algorithm must be terminated, as there is no strategy for following iterations. This drawback further restricts the choice of initial guesses to a set of parameters with which the forward model can be solved.

Computational cost
The computational cost of the algorithm is dominated by the cost for the forward model evaluations. Therefore it scales at least linearly with increasing number of parameters n x , as more simulations are required per iteration for the finite difference approximation of the Jacobian. Additionally, a higher dimension in parameters of the inverse problem leads to more complex objective functions and therefore potentially longer search paths. In the same sense a higher number of parameters n x will potentially lead to a smaller region around the optimum for the individual parameters from which the initial guess has to be chosen to be able to find an optimum.
In presented examples the algorithm converged or failed in 20-40 iterations.

High number of parameters
From all of the above points it becomes clear that the presented method does not scale arbitrarily well for high number of parameters. The more parameters involved in the model, the more likely it is, that they are differently significant to the forward model solutions and possibly interact in their contribution to the objective function. The objective function gets more complex with higher number of parameters and the probability, that it becomes multimodal, i.e. more than one local optimum exists, increases. The more complex the objective function is, the more likely it is, that the search path leads to parameter regions, where the forward model cannot be evaluated and the inverse analysis yields no result. Higher number of parameters also leads to higher computational cost of a finite difference approximation and search paths grow longer in more complex objective functions increasing the number of finite difference approximations necessary. Overall it is expected that the cost and general applicability of the proposed method scale poorly for high parameter dimensions. For example the inverse analysis of subdomain shapes and therefore an element-wise definition of material parameters can be considered as very high dimensional in general. To overcome the described problems with high dimensions so called Patched Basis Functions can be defined to find the subdomain patches [44].

Importance of initial guess
The presented algorithm includes a local optimizer and it cannot find a reasonable solution for every given combination of parameters as initial guesses. The availability of a good initial guess is crucial, as it decides if the method converges and also how many mostly costly forward model evaluations are required. Nevertheless the presented algorithm itself can help to find a good initial guess, when it is used with a reduced set of parameters. In the application with real experiments it is unknown if an optimum that was found with one initial guess using the presented approach is a global optimum with regard to given data. Hence inverse analysis results must always be carefully interpreted with respect to plausibility of the results. In doubt it is always possible to validate results by using a second and significantly different initial guess and compare the results.

Model selection
Deterministic inverse analysis is used to find a point estimate only, representing a local best fit of parameters in a chosen forward model. It does not provide a general answer to what type of measurement error is inherent in the data and also not if the forward model used for optimization is itself a good choice or should be improved. This statement is not restricted to the material model but also holds for the physical model and the boundary conditions used for the forward model. In the case of raw flow cell experiment results for example, the best choice for a specific hyperelastic material model is unknown and a Saint-Venant-Kirchhoff material model with its two parameters is just one very simple choice. The selection of a suitable forward model is up to the analyst with presented approach, but also could be included in the inverse analysis by identifying the best parameters for different models and comparing the overall approximative error (like also done in e.g. [34]). If one does that it is very important to relate the approximative quality to the number of parameters in the model (e.g. by the Bayesian or the Akaike information criterion) as we are seeking a predictive model and not just a fit to some data points (see e.g. [34] or [45]).

OCT imaging
The presented measure used for the objective function works independent of the physical model and the spatial dimension of the model. An important feature is that it can easily be used tailored to the data gained from OCT. For example, in the unloaded state a three dimensional scan of a biofilm can be used to build a mesh for computational evaluation and the objective function can be defined with so called B-Scans from the loaded and therefore deformed biofilm, whenever speed of image acquisition is the limit to measurement quality. B-Scans are two dimensional plane scans of the flow cell [7] and naturally faster to acquire than three dimensional images which are just stacks of multiple B-Scans.

Combination of data
Only examples that include information of the initial non deformed geometry of the biofilm incorporated in the forward model mesh and one quasi steady state solution are shown for the sake of concise presentation. But obviously also applications with model evaluations in different time steps or load steps to asses nonlinear material parameters in more complex material models or viscous effects of the biofilm can be easily treated. In such applications, questions of scales and weights of the contribution to the residual r in the objective function must be answered, to not blindly weigh highest displacements the highest. Nevertheless it has been shown in [36] that also results from different experiments on the same specimen can be scaled and used in a single Levenberg-Marquardt based inverse analysis.

Time scales
In the application of the method, the experiments and results should be grouped in meaningful sets by the time scales involved. One use case can be to first analyze a deformation experiment in the flow cell and determine material parameters in a simple material model and as a second step use other experiments with a more extensive model like additional growth and determine the growth parameters (like in the growth example). In this scenario the growth and hyperelastic material parameters can be regarded as decoupled as time scale for growth is orders of magnitude larger than the time scale for fluid-solid interaction under constant fluid inflow rate.

Convergence criterion
A convergence criterion based on combined maximal number of algorithm iterations and a maximum gradient based error presented itself as a good choice. Even with the presented artificially generated data it was not easily predictable how low the remaining residual error between measured deviation of forward model outcome and experiment observation is. Especially if artificial noise was added to the data, it showed that there was an individual distinct level of residual error, when no better set of parameters could be found. This is also obvious as the noisy data does not define a real physical solution and hence an error has to show up. This combination of convergence criteria can therefore be recommended and needs to be tuned to the specific application depending on the unknown structure of uncertainties in the experiment result data and the cost of forward model evaluations.

Conclusion
An inverse analysis method for biofilms has been presented and successfully tested on several selections of significant parameters in different aspects and models of biofilm physics. The algorithm works with a local best fit for parameter estimates in given forward models related to experimental results in the sense of least squares. A simple hands-on measure for the comparison of shapes of biofilms in deformation experiments has been presented and tested. It has been shown for a variety of different meaningful models for biofilms, like homogeneous, heterogeneous, poroelastic and biofilm models including surface growth, that the presented approach allows the inverse analysis for multiple key parameters at once. As the presented measure for the difference of a biofilm surface between experiment and model evaluations is based purely on their shapes, the approach can without restrictions be used for further different physical models, like ones for detachment, self contact and viscoelasticity.
Author contributions HW implemented the inverse analysis approach, run the simulations, prepared the figures and primarily drafted the manuscript. WAW worked out the general conception of the project. Both authors contributed to the discussion of results and prepared and approved the final manuscript.

Funding
Open Access funding enabled and organized by Projekt DEAL. German Research Foundation (DFG) with project number WA 1521/22.

Availability of data and materials
The research code is hosted on a private GitLab repository at Leibniz Rechenzentrum (LRZ) in Garching. The generated numerical results and digital data are held on machines that are backed up on servers managed by the LRZ in Garching. The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.