Recovering the Fragmentation Rate in the Growth-Fragmentation Equation

We consider the inverse problem of determining the fragmentation rate from noisy measurements in the growth-fragmentation equation. We use Fourier transform theory on locally compact groups to treat this problem for general fragmentation probabilities. We develop a regularization method based on spectral filtering, which allows us to deal with the inverse problem in weighted ${L}^2$ spaces. %Our approach regularizes the signal generated by differential operators in the frequency domain. As a result, we obtain a regularization method with error of order $O(\varepsilon^{\frac{2m}{2m+1}})$, where $\varepsilon$ is the noise level and $m>0$ is the {\em a priori} regularity order of the fragmentation rate.


Introduction
Growth-fragmentation equations describe, in a quantitative way, the evolution process of the density of an ensemble of particles that divide according to a certain fragmentation rate.In other words, we assume that each particle grows over time, and splits into smaller particles, in such a way that the total mass is conserved.This model is used in biological phenomena for instance, in cell division processes [1,2,3,4] and protein polymerization [5].It also appears in telecommunications [6].The goal of this article is to address the inverse problem for a class of such equations making use of the asymptotic behavior of their solutions which in turn, after proper scaling, converge to a stable distribution in time [1,7,8].
The growth-fragmentation equation considered here can be expressed using the following integro-differential equation, where n(t, x) represents the density of particles of mass x at time t, with initial condition n 0 (x).The function g(x) is the growth rate for particles of mass x.The function K(x|y) represents the conditional probability that a particle of mass y splits into k smaller particles of mass x.The function B(x) is the fragmentation rate of particles of mass x.
A natural question concerns to the asymptotic behaviour of the population density n(t, x), when t → ∞.In [7,8], it is shown that under certain conditions on the coefficients g, B, and K, there exists an unique eigenpair (N, λ) such that ( Moreover, in appropriate weighted norms, we have that n(t, x)e −λt → N(x), when the time t goes to infinity.
Thus, we shall focus on the pair (N(x), λ), which corresponds to the eigenpair of a linear operator [7].They have the biophysical interpretation of the asymptotic distribution of the population density and the corresponding convergence rate.In principle, such quantities could be experimentally observed in some situations due to the exponential decay rate.This was used for instance in [9].The function N alluded to above is the stable distribution which is observed in a broad class of structured population models [1].
The importance of studying N(x) instead of n(t, x) is that the variable t is removed, which reduces the study of a two-dimensional problem to a one-dimensional.
We are interested in how to recover in a stable way the fragmentation rate B(x) from noisy measurements of N(x) and λ.Our strategy is to study H = BN rather than B, as proposed in [10,11], and then use truncated division by N to recover B. Thus, the inverse problem under consideration is how to recover in a stable way H from approximate knowledge of N and λ, in which H is the solution of where In [10,12,9] the inverse problem was treated for the equal mitosis case, that is, when a particle splits into two small particles with half of its mass.In this case, the parameters are A more general class of conditional probability density kernels is given by the self-similar ones, that is, when where P is a probability measure in [0, 1].Self-similar probabilities arise when the fragmentation depends on the ratio between the size of the mother particle and the size of the next-generation particle.The inverse problem for self-similar probabilities was treated in the works [11,13].
The aim of the present work is to address this inverse problem for more general probabilities.Namely, we consider conditional probabilities of the form where ρ is an increasing diffeomorphism from R to R + , and P is a probability on (0, ρ(0)).
The transport operation ⊙ ρ is defined by We shall denote conditional densities in (7) as transport probabilities.Observe that if we take the exponential function ρ(x) = e x , then the transport operation for this function is the usual product on the positive real line.Hence, transport probabilities generalize self-similar probabilities.
We treat the inverse problem described in Eq. ( 3).Here, we do not require that the mass must be conserved in the division process.This is in order not to restrict this method only to biological models, and also to use it in real-world applications.treating this problem from a more general point of view.We do that in two steps.Firstly, we guarantee that under some assumptions in N, there exists a unique H solution of Eq. (3).
To do that, we establish that on proper spaces the operator kK − Id is an isomorphism.
The second step is to guarantee the stability of the inverse problem.For that, we propose a new regularization method.This method is based on the implementation of the Fourier transform theory, and spectral filtering techniques.
We show that the Quasi-Reversibility approach [12] is a particular case of our method.
Compared to regularization methods based on convolution [11,13], our method improves the error order.Namely, we obtain an error of order O(ε ), where ε is the noise level.This article is organized as follows: in Section 2, we study the transport operation and some relations with the Fourier transform on locally compact abelian groups.In Section 3, we discuss the invertibility of the operator kK − Id in proper spaces.In Section 4, we present a new regularization method to treat the stability of the inverse problem.In Section 5, we give examples for which some of the Hypotheses Eq. (3.1) or Eq.(4.1) are satisfied.Finally, in Section 6, we present the numerical implementation of our method.

Preliminaries
To fix the notation, we briefly review some concepts of Fourier transform for locally compact groups.We refer the reader to [14] for more details.We recall that the Haar measure on an abelian group G is the unique non-negative and regular measure, up to a positive multiplicative constant, which is translation invariant.Here, by translation, we mean the action of multiplication by group elements.
We define the transport operation ⊙ ρ as where ρ is an increasing diffeomorphism from R to R + .The pair (R + , ⊙ ρ ), equipped with the induced topology of R, is a locally compact group.Moreover, the Haar measure µ ρ is given by the measure where dx is the Lebesgue measure on the positive real line.In this context, we can develop the theory of the Fourier transform on the group (R + , ⊙ ρ , µ ρ ).
The Fourier transform F ρ on the group (R + , ⊙ ρ ) is defined for a function f (x) ∈ L 1 (R + , µ ρ dx) and a real number ξ as In this general context, the Fourier transform theory on (R + , ⊙ ρ ) can be developed in a similar fashion as in the standard case (R, +).In particular, we obtain the inversion and the Plancherel theorems.See [14] for more details.
Theorem (Inversion theorem [14]).Suppose that f ∈ L 1 (R + , µ ρ dx) and then for a.e positive number x we have Theorem (Plancherel Theorem [14]).The Fourier transform F ρ , restricted to ) is an isometry (with respect to the L 2 norm) onto a dense linear subspace of L 2 (R, dx).Hence, it may be extended, in a unique manner, to an isometry from L 2 (R + , µ ρ dx) to L 2 (R, dx).
In order to extend the Plancherel theorem, we define the L 2 ρ ,s space.They are a natural generalization of the L 2 (R + , µ ρ dx) space.More precisely, the L 2 ρ ,s space is defined as L 2 (R + , e 2πsρ −1 (•) µ ρ ).Observe that there exists a canonical isomorphism T s between L 2 ρ ,s and L 2 ρ ,0 = L 2 (R + , µ ρ dx), which is given by We now define the Fourier transform F ρ ,s on L 2 ρ ,s as Thus, by the Plancherel Theorem, the generalized Fourier transform F ρ ,s is an isometry.
For the case when P is a probability measure on R + satisfying we define the Fourier transform F ρ ,s of the probability P by for any real number ξ.

The Fourier transform and the convolution operator
We now consider the convolution-type operator induced by the probability P and given by In the same manner, as in the case of R, and under some conditions of integrability, the Fourier transform F ρ ,s of a convolution operator becomes a multiplicative operator.
Property.Let P be a probability measure satisfying Equation (15), then the convolution operator defined in Eq. ( 17) satisfies the multiplication property for all f in L 2 ρ ,s .

The Fourier transform of d dx gN.
The Fourier transform on the real line has the property of diagonalizing differential operators.We now state an analogous version of this fact for the Fourier transform F ρ ,s .
Proposition 2.1.Suppose that S belongs to the Sobolev space H 1 0 (0, ∞), and that the functions are in L 2 ρ ,s .Then, we have Proof.Using integration by parts, and the smoothness assumption of the function S, we have that and this proves our result.
3. Invertibility of the operator kK − Id.
Let K : L 2 ρ ,s → L 2 ρ ,s be the operator defined as We use the ideas of [11] in order to guarantee the invertibility of kK − Id.We prove that under additional assumptions on the probability P , the operator kK − Id, has a bounded inverse in some subspaces of L 2 ρ ,s .If we apply the Fourier transform and the Eq. ( 18), we get that kK − Id is in fact a multiplication operator Thus, the operator kK − Id has a bounded inverse on L 2 ρ ,s if the function is bounded from below by positive constant on the real line.In some cases, the function |k F ρ ,s P (•) − 1| never vanishes, but goes to zero when x converges to 0 or ∞, in this situation, the function |k F ρ ,s P (•) − 1| is bounded from below on compact sets of the real line.Moreover, for computational purposes, we are interested in a reconstruction on compact intervals.In this case, we focus on the invertibility of the operator on compact sets.To find a bounded inverse on an open set U, we assume that the probability P satisfies the following hypothesis.In Section 5 we give examples where the previous hypothesis is satisfied.
We consider the Paley-Wiener spaces P W ρ ,s (U), as the subspace of L 2 ρ ,s of all the functions f whose Fourier transform F ρ ,s has support on U. Observe that if f is on P W ρ ,s (U), then Kf also lies on P W ρ ,s (U).Thus, the operator kK − Id from P W α ρ ,s to itself is well defined.Using the Fourier transform F ρ ,s , and Eq. ( 18), we conclude the following proposition, which guarantees the existence of a bounded inverse.Proposition 3.1.Suppose that the probability P satisfies Hypothesis 3.1, then the operator kK − Id from P W ρ ,s (U) to itself has a bounded inverse.

Regularization of the inverse problem
In this section, we discuss some methods to regularize the inverse problem, i.e, recover H = BN from noisy measurements of N in some L 2 ρ ,s norm, in such a way that we can control the approximation error.
The main difficulty arises from the fact that we cannot estimate the L 2 ρ ,s norm of d dx gN from N .The principal idea in this section is to use the Fourier transform to turn the differential operator d dx into a Fourier multiplication operator [15], and then, we regularize the inverse problem using spectral filters as described in [16].

Regularization by spectral filtering
We now present our strategy to regularize the inverse problem as follows.Using the Proposition 2.1, we see that under the Fourier transform, the differential operator is simply a quasi-multiplication operator, whose multiplication part is given by (2iξ − s)π.
Unfortunately, this function is unbounded in non-compact sets.Then, noisy data with small errors can approximate solutions with large errors, that is, it has a de-regularizing effect.To regularize the inverse problem, we define a multiplication approximation for d dx gN, for which the multiplier is given by (2iξ − s)π by a regularized filter f (α, ξ).To be more precise, we define the approximation H α (N ε , λ ε ) by where and let H = H(N, λ) be the unique solution of Eq. ( 3) in P W ρ ,s (U).The following theorem Preprint 9 gives error estimates for H α (N ε , λ ε )−H(N, λ) .It was inspired by [10], yet, the method is different because it focuses on regularizing the spectral signal of the differential operators.
Theorem 4.1 (Spectral regularization).Let K be a subspace of P W ρ ,s (U).Suppose that there exists a positive constant C, such that for all α ∈ [0, C) the bilinear operator where the constant M does not depend on α.Then, we have the estimate Proof.By the triangle inequality, we obtain which implies the result.
To apply the above result, we need to guarantee that for a fixed α, the operator H α is well defined and bounded.For simplicity, we first consider the case when and are bounded functions on U. To show that the operator H α is well defined, we consider the subspace D(U) of P W ρ ,s (U), which consists of all functions N such that are in P W ρ ,s (U).In that case, it is straightforward to show that the operator H α from U) is well defined, and bounded.We use modified versions of Tikhonov and Landweber filters, to show that all the functions in Equation (25) are bounded.These filters are commonly used in regularization theory for compact operators [16].

Tikhonov filtering
In the classical case of linear operators, the solution x to the problem Ax = y, where A : H → H is a non-negative linear and compact operator, can be regularized by using the filter in the spectral variable.Indeed, considering the singular value decomposition of A given by the regularized solution x α to the problem Ax = y is given by See [16] for more details.For the problem under consideration, we modified the above function and consider the filter We assume that the probability P satisfies Hypothesis 3.1.Using these filters, we observe that the functions in Equation ( 25) are bounded with respect to the variable ξ.Under these assumptions, the bilinear operator H α is bounded.In fact, a straightforward computation shows that for some positive constant C, which only depends on s, P , g, ρ −1 .
We now state the following regularization method, based on the filter functions defined above.
Then, we can choose the optimal parameter α = √ ε to conclude Proof.Using the Fourier transform in Eq. ( 3), together with Proposition 2.1, we have Therefore, we obtain for some constant M ′ .Hence, by Theorem 4.1 we obtain the result.

The quasi-reversibility Tikhonov filtering
The quasi-reversibility method proposed in [12] regularizes the inverse problem for the case of equal-mitosis.The idea of quasi-reversibility goes back to the work [17].This method is based on adding a perturbation of a differential operator.In [13], an extension of this method was proposed for more general probabilities.We now present a different generalization using spectral filters.In fact, we show that there exists a relation between the quasi-reversibility method and our approach in the case of self-similar probabilities.
Let us first describe our method without going into technical details.For each real number j, we defined the following modified Tikhonov filters If the quotient factor blows up, then the last equation is not well defined.To avoid this kind of problem, we assume that the probability P satisfies a stronger hypothesis than 3.1, namely, we assume that In Section 5, we give some examples for which the above hypothesis is satisfied.If we assume that N is a smooth function satisfying Condition 29, we see that under the Fourier transform F ρ ,s , the approximation in Eq. ( 23) solves the following differential equation where Thus, our method can be seen as a perturbation method [18].For the self-similar case, that is, when ρ(x) = e x , we see that Taking j = 1, the perturbation method defined in Eq. ( 27) reduces to the quasi-reversibility method proposed in [13].Using the same ideas of the proof of Theorem 4.2, we prove the next result.
for some constant M which depends on P , g, ρ −1 , and N.

The Landweber's method
Using the Tikhonov and quasi-reversibility filters, we obtain an approximation error of order O( √ ε).We improve this error order, and thus we obtain a better approximation.
For that, we use a new filter, which is a modification of the classical Landweber filter To implement these filters in Theorem 4.1, we first establish the following inequalities.
Lemma 4.1.There exists a positive constant A, such that for all x ∈ R >0 and for all positive α, the following estimate holds where A depends on s, P , g, ρ −1 .Moreover, for all m ∈ R, such that m < 2α, the following inequality holds Proof.Let us prove the first inequality.Using the change of variables z = 1 − 1 1+x 2 , we see that Since z ∈ [0, 1], then applying the mean value inequality to the function z α , we obtain To prove the second inequality, we observe that at the point ω = 2α−m 2α , the function (1 − z) m z 2α−m attains its maximum for x ∈ [0, 1].Thus, for all x ∈ R we have As a consequence of the previous lemma, we have that the functions in Eq. ( 25) are bounded.Thus, the bilinear operator H α is bounded.In fact, a straightforward computation shows that for some positive constant A, which only depends on s, P , g, ρ −1 .Now, we apply Theorem 4.1 to the Landweber filter [16].For that, we require a smoothness condition of order m > 0 for the function N Then, for all noisy measurement (N ε , λ ε ) of (N, λ) in D(U) × L ∞ (0, ∞), the approximation H α (N ε , λ ε ) defined in Eq. ( 23) using the spectral filter of Eq. (28) satisfies the error estimate for some constant M which depends on s, m, P , g, ρ −1 , and N. Suppose that (N ε , λ ε ) √ m.We can choose the optimal parameter given by , which satisfies m < 2α, to conclude Proof.Using the Fourier transform in Equation (3), and Proposition 2.1, we have that Next, we apply the second estimate of Lemma 4.1 to obtain hence by Theorem 4.1, we obtain the result.

The unbounded case
In many cases, the functions in Equation ( 24) are not bounded.For instance, in the selfsimilar case.We now extend our results to the unbounded case.The idea is to regularize the functions in Equation (24) using a spectral filter.To do that, we write and , where α ∈ R >0 , is the regularization parameter of the functions (24).Observe that the exponential decay guarantees that using g α in the place of g, then the functions in Eq. ( 24) are bounded.Now, we show that if the function N has fast decay, then the bounded function g α can be used to regularize the inverse problem, even if g is not bounded.
Proposition 4.1.Assume that the function N satisfies We define H α as the solution of Eq. (3), where the solution is associated with the function g α instead of g.Then, we have that where C is a constant, which depends on the operator K and the number k.
Proof.Observe that then, using the Fourier transform F ρ ,s , we obtain Thus, the function H α is a controlled approximation for H.In this case, the function H α is the solution of Eq. ( 3), associated with the bounded function g α .Therefore, we can use some of the previous methods (Tikhonov or Landweber) to regularize H α , and by Proposition 4.1 to regularize H.

Examples
To use the previous regularization methods, we need to verify that the probability P satisfies Eq. ( 15), and some of the Hypotheses 3.1 or 4.1.In this section, we study some examples which satisfy the previous conditions.First, we discuss examples arising from the self-similar probabilities, that is, when ρ(x) = e x and k = 2.
To check Hypothesis 4.1, it is sufficient to consider the case j = 0.In fact, if Hypothesis 4.1 holds for j = 0, then by triangle inequality we see that for all j = 0 Thus, the Hypothesis 4.1 holds for all real numbers j.We now study for which values of s, and open sets U, the Hypothesis 4.1 holds.Without loss of generality, we assume that j = 0.

The equal-mitosis
An important example of a self-similar case is the equal-mitosis.In equal-mitosis the conditional probability is given by P = δ x= 1

2
. In this case, the Fourier transform of P is given by .
Thus, by the triangle inequality Therefore, we conclude that this probability satisfies the Condition (15) and Hypothesis 3.1, for all s ∈ R \ { 2 π }.In order to verify Hypothesis 4.1, we consider the following cases.
For such s, the Hypothesis 4.1 is satisfied for all ξ ∈ R. In fact, by the triangle inequality • Second case: s < 2/π.
In this case, the Hypothesis 4.1 holds for each open bounded set, and α small enough.
In fact, using triangle inequality we have Then, using estimate 30, we conclude that for α small the expression |F ρ ,s (P )(ξ) − 1| − 2πiαξ is positive.
Thus, we can apply our methodology to deal with the inverse problem for these cases.In Section 6, we will develop some numerical simulations for the previous cases and see the effectiveness of our approach.

Numerical solution of the inverse problem
In this section, we recover numerically the solution of the inverse problem using the regularization methods proposed in Section 4. That is, we recover H from noisy measurements of N. If we assume that the noisy measurement N ε is a smooth function satisfying Condition 29, then there exists a unique solution H ε α for the Equation (3) associated with N ε .The purpose of this section is to explore numerically how close is H ε α to H, when the noisy data N ε is close to N. To do that, we first construct an approximation for N, using the numerical schemes proposed in [10,13], and then, we add random noise to the data N.We construct the approximation H ε α using Equation (23), where α is the optimal regularization parameter for each method.

Parametric specification.
We now present some numerical simulations for the equal-mitosis case.Here, the parameters are s = 0, k = 2 and ρ(x) = e x , and the probability distribution is given by To guarantee the existence and uniqueness of the solution of the direct problem, and also the condition Eq. (24), we select fragmentation and growth rates with fast decay at 0 and ∞.See [19].To be more specific, we use the growth rate g(x) = xe −(x+1/x) and the fragmentation rate In this experiment, we use the space L 2 ρ,0 to recover the function H on (0, 10].We also use the parameter j = 1 for the quasi-reversibility method and m = 10 for the Landweber filter.

Construction of N.
We construct the function N on the interval (0, 3] using the numerical scheme proposed in [10,13].Here, the initial condition is given by We use a regular grid on (0, 3] with 500 points.For the evolution process, we use a regular grid on (0, 200], with 10 4 points.We plot the function N in Figure 1.

Reconstruction of H and B.
We now consider a noisy approximation (N ε , λ ε ) for the eigenpair (N, λ), obtained by adding a random noise to the data.That is, we assume that where R ε , S ε are random noises uniformly distributed in [−ε, ε].We recover the approximation H ε α using the noisy measurement (N ε , λ ε ).We plot the approximation H ε α for different noise levels ε.Here, we use the values ε = 10 −2 (Figure 2), and ε = 10 −3 (Figure 3).The parameters used are m = 10 for the Tikhonov method and j = 1 for the Landweber method.
To recover the fragmentation rate B from H = BN we use the truncate division by We compare the numerical error of the reconstructions of the functions B and H for small values of ε.Here, we use the L 2 [0, 3] norm to estimate this error.This norm is computed using rectangular integration.For this experiment, we assume that ǫ ∈ [10 −4 , 0.5].We plot the error in logarithmic scale in Figure 4. Observe that for small values of ǫ the Landweber filter gives a better reconstruction, which is following our theoretical results.

Conclusion
In this article, we treated the inverse problem for the growth-fragmentation equation of a wide class of transition probabilities.Namely, we deal with transport probabilities which are generalizations of self-similar probabilities.We developed a new approach to regularize the inverse problem associated with the transport probabilities.Our approach is based on the Fourier transform theory for locally compact groups.We regularized the Fourier transform of differential operators using several filters in the spectral variable, such as modifications of the Landweber and Tikhonov filters.
For each method, we obtained their respective error estimates.The Landweber method provides an algorithm to recover the fragmentation rate B, with order O(ε 2m 2m+1 ), where m is the degree of smoothness of B, as proved in Theorem 4.4.The error obtained using the Landweber method is better compared with other methods.This fact was verified by numerical simulations.
Our theoretical approach has been focused on transport operation induced by diffeomorphisms.A natural continuation of this work is to deal with transport operations induced by arbitrary functions.Another possible follow-up is to apply our methodology to problems arising in data science, as well as biological problems as in [20,21].
One potential application of the problem under consideration concerns modeling normal prion protein (PrP(C)) and infectious prion protein (PrP(Sc)) populations interacting in an infected host [5].Indeed, Perthame and Doumic in [22] proved key asymptotic results for the stable population where our results could be applied.
We note that even in the particular zero-growth case, g ≡ 0, we have an interesting topic for further exploration since it allows for some simplifications in the operator of Equation ( 3).In such cases Doumic et al. [23] have addressed the issue of reconstructing the conditional probability kernel under a self-similar assumption and knowledge of the fragmentation rate B. This shall be addressed in a future publication.

Hypothesis 3 . 1 .
There exists an open set U, such that the function |k F ρ ,s P (•) − 1| is bounded from below by a positive constant on U.

Theorem 4 . 2 (
Tikhonov regularization).Assume that the probability P satisfies Condition 15, and Hypothesis 3.1.Moreover, we assume that N ∈ D(U) satisfies all the assumptions of Proposition 2.1, and

Theorem 4 . 3 (
Quasi-reversibility method).Assume that the probability P satisfies Eq. (15) and Hypothesis 4.1.Moreover, assume that N ∈ D(U) satisfies all the assumptions of Proposition 2.1.If the functions ξ F ρ ,s ( d dx gN) and ξ F ρ ,s ( d dx N), are in L 2 (R), then, for all noisy measurement

Figure 1 :
Figure 1: Construction of N, solution of the direct problem for the equal-mitosis case.

Figure 4 :
Figure 4: Numerical error of the reconstruction of the functions B and H using the L 2 norm for several values of ǫ.Note that errors are on a logarithmic scale