Structured Models for Cell Populations : Direct and Inverse Problems

Structured population models in biology lead to integro-differential equations that describe the evolution in time of the population density taking into account a given feature such as the age, the size, or the volume. These models possess interesting analytic properties and have been used extensively in a number of areas. In this article, we consider a size-structured model for cell division and revisit the question of determining the division (birth) rate from the measured stable size distribution of the population. We formulate such question as an inverse problem for an integro-differential equation posed on the half line. The focus is to compare from a computational view point the performance of different algorithms. Finally, we will discuss real data reconstructions with E. coli data. 1 Structured Population Models Let us consider a population density n = n(t, x) where t ∈ R+ is time and x ∈ R+ is a parameter, e.g., the mass of a unicellular organism, the DNA content of a cell, the cell age, or its volume. A more specific example is the cell volume in a population of E. coli. This unicellular organism is commonly found in the lower intestine of warm-blooded organisms. The cells are typically rod-shaped and are about 2μm long, with 0.5μm diameter. Then, the average volume falls between 0.6 − 0.7μm3. These unicellular organisms are extensively studied in vitro and in vivo. Let n(t, x) denote the population density of cells of “size" x at time t. Then, n satisfies the integrodifferential equation ∂tn(t, x) + ∂x[g(x)n(t, x)] = ∫ ∞ 0 k(x, x′)n(t, x′)dx′, (1) where g(x) is the microscopic growth rate of individuals at size x, and k(x, x′) is the proportion of cells of size x′ that divide into two cells, one of size x, and the other of size x′ − x. Under this generality, the model is hard to calibrate and to make predictions. In one spatial dimension, the amount of data needed to estimate the division kernel k(·, ··) would be unrealistic in practical applications without further hypothesis. See [1–3] for variants of the division ae-mail: vvla@impa.br be-mail: zubelli@impa.br DOI: 10.1051/ C © Owned by the authors, published by EDP Sciences, 2015 / 0 0 ( 2015) 201 conf Web of Conferences ITM itm , 000 5 0 05 5 16 16 This is an Open Access article distributed under the terms of the Creative Commons Attribution License 4.0, which permits distribution, and reproduction in any medium, provided the original work is properly cited. 2 Article available at http://www.itm-conferences.org or http://dx.doi.org/10.1051/itmconf/20150500016 kernel estimation problem. Thus, in what follows, we shall consider the following simplified version of the model given by Equation (1): { ∂tn(t, x) + ∂xn(t, x) + B(x)n(t, x) = 4B(2x)n(t, 2x), x, t ≥ 0, n(t, x = 0) = 0, t > 0, n(0, x) = n0(x) ≥ 0, x ≥ 0. (2) We remark that the choice of g ≡ 1 was made and that a natural alternative would be that of an affine function. It is well-known [4–6] that there exist unique λ0 and N = N(x), forming an eigenpair, such that, after a suitable time re-normalization, the solutions of (2) satisfy the limit n(t, x)e−λ0t −→ ρN(x), as t → ∞, (3) where the limit is considered under weighted Lp topologies, and the pair (λ0,N) is a solution for { ∂xN(x) + (λ0 + B(x))N(x) = 4B(2x)N(2x), x ≥ 0, N(x = 0) = 0, N(x) > 0, for x > 0, ∫ ∞ 0 N(x)dx = 1. (4) Moreover, exponential rates have been proved to hold for fairly general rates B. See [6, 7]. Such N is the so-called stable-size distribution. Let us assume that the birth rate B is a measurable function and such that, there exist positive constants Bm, and BM satisfying 0 < Bm ≤ B(x) ≤ BM < ∞. (5) Then, we can define the direct problem as, given a birth rate B satisfying such conditions, finding the eigenpair (λ0,N) of Problem (4). The corresponding inverse problem is to recover in a stable way the birth rate B from noisy data N and the rate λ0. The first step is to study the so-called parameter-to-solution map. It can be shown that the following result holds. Theorem 1 (Perthame-Zubelli [8]) Under assumption (5), the map B → (λ0,N), from L∞(R+) into [Bm, BM]× L1 ∩ L∞(R+) is continuous under the weak-∗ topology of L∞(R+), locally Lipschitz continuous under the strong topology of L(R+) into L(R+), and of class C1 in L(R+). 2 The Inverse Problem If the measurement N were very smooth, one could directly consider to solve for B, the cell-division equation (4) written with y = 2x, 4B(y)N(y) = B (y/2) N (y/2) + λ0N (y/2) + 2∂yN (y/2) , y > 0. (6) This is a well-posed equation on B as long as N satisfies regularity properties such as ∂yN (y/2) is in Lp, for some p ≥ 1. However, this is not the case for reasonable practical data. One possible regularization strategy is the following α∂y(BαN) + 4Bα(y)N(y) = Bα (y/2) N (y/2) + λ0N (y/2) + 2∂yN (y/2) , y > 0, (7) with BαN(0) = 0, where 0 < α < 1 is a small parameter. In [8], it was shown that Equation (7), when solved for B has desired properties such as strong stability with respect to perturbations. As a consequence, given a noise level ε, choosing α = α(ε), we have ‖Bε,α − B‖L2(N2 εdx) ≤ K √ ε. ITM Web of Conferences


Structured Population Models
Let us consider a population density n = n(t, x) where t ∈ R + is time and x ∈ R + is a parameter, e.g., the mass of a unicellular organism, the DNA content of a cell, the cell age, or its volume.A more specific example is the cell volume in a population of E. coli.This unicellular organism is commonly found in the lower intestine of warm-blooded organisms.The cells are typically rod-shaped and are about 2μm long, with 0.5μm diameter.Then, the average volume falls between 0.6 − 0.7μm 3 .These unicellular organisms are extensively studied in vitro and in vivo.
Let n(t, x) denote the population density of cells of "size" x at time t.Then, n satisfies the integrodifferential equation where g(x) is the microscopic growth rate of individuals at size x, and k(x, x ) is the proportion of cells of size x that divide into two cells, one of size x, and the other of size x − x.Under this generality, the model is hard to calibrate and to make predictions.In one spatial dimension, the amount of data needed to estimate the division kernel k(•, ••) would be unrealistic in practical applications without further hypothesis.See [1][2][3] for variants of the division a e-mail: vvla@impa.brb e-mail: zubelli@impa.br We remark that the choice of g ≡ 1 was made and that a natural alternative would be that of an affine function.
It is well-known [4][5][6] that there exist unique λ 0 and N = N(x), forming an eigenpair, such that, after a suitable time re-normalization, the solutions of (2) satisfy the limit where the limit is considered under weighted L p topologies, and the pair (λ 0 , N) is a solution for Moreover, exponential rates have been proved to hold for fairly general rates B. See [6,7].Such N is the so-called stable-size distribution.
Let us assume that the birth rate B is a measurable function and such that, there exist positive constants B m , and Then, we can define the direct problem as, given a birth rate B satisfying such conditions, finding the eigenpair (λ 0 , N) of Problem (4).
The corresponding inverse problem is to recover in a stable way the birth rate B from noisy data N and the rate λ 0 .The first step is to study the so-called parameter-to-solution map.It can be shown that the following result holds.

The Inverse Problem
If the measurement N were very smooth, one could directly consider to solve for B, the cell-division equation ( 4) written with y = 2x, This is a well-posed equation on B as long as N satisfies regularity properties such as ∂ y N (y/2) is in L p , for some p ≥ 1.However, this is not the case for reasonable practical data.
One possible regularization strategy is the following with B α N(0) = 0, where 0 < α < 1 is a small parameter.
In [8], it was shown that Equation (7), when solved for B has desired properties such as strong stability with respect to perturbations.As a consequence, given a noise level ε, choosing α = α(ε), ITM Web of Conferences In this note we take a more pedestrian approach to tackle the inverse problem than that used in [8][9][10].More specificaly, we want to find minimizers for the following Tikhonov-type functional with B ∈ L 2 (R + ), satisfying (5), and α = 0.05.The penalization functional used are: , and Kullback-Leibler: where with Kullback-Leibler, B 0 (x) = 0.01 with E. coli data and B(x) = 1 with synthetic data, and with the smoothing functional, B 0 equals to the division-rate obtained in [9] with E. coli data and also B(x) = 1 with synthetic data.

Numerical Results
We now present some numerical examples of the calibration of the division-rate B from the stable distribution N using Tikhonov-type regularization and statistical inverse-problem techniques, with synthetic as well as real E. coli data.The latter data is the same used in [9].The numerical solution of the PDE problem ( 4) is obtained with the numerical scheme presented in [10].
By describing the direct problem under a discrete Bayesian model [11], we build a posterior density, which is proportional to the exponential of the functional in Equation ( 8).The reconstructions are obtained with Maximum a Posteriori (MAP) and Conditional Mean (CM) point estimators.The MAP estimate as well as the Tikhonov-type reconstructions are obtained with the MATLAB's function LSQNONLIN.To evaluate the CM, we use the Metropolis-Hastings algorithm, which is initialized with the MAP estimate and was generated with 10 4 samples.See [11].With synthetic data, the minimization starts with B(x) = 1, and with E. coli data it starts with the division rates obtained in [9].
The synthetic data is generated with smooth and non-smooth division-rates B, in a mesh two times finer than the one used to solve the inverse problem.The data is corrupted by a multiplicative Gaussian noise with mean one and covariance matrix δ 2 Id, where δ = 0.05.Figures 1 and 2 present the reconstructions of B with synthetic data.In both cases, reconstructions are compatible with the true B. Note that, for the non-smooth B, we also used a l 1 prior density on the first derivative of B in the statistical approach.For the smooth case, instead, we used the l 2 prior density on the second derivative of B.
Figures 3 presents the reconstructions of the division-rates using Tikhonov regularization and the statistical approach interpolated by cubic spline, both compared with B obtained in [9]. Figure 4 Hybrid and Multiscale Modelling in Biology 00016-p.3 presents the corresponding stable density functions for both cases.Note that, we used the E. coli data without interpolation.Visually, our reconstructions present a similar shape as the one obtained in [9].Moreover, the corresponding stable size-distributions presented a good adherence to the experimental data.

Discusion
The solution of the inverse problem relies heavily on the existence of a limiting behavior for population structure.Indeed, by looking at the long time behavior, we do not need to assume a particular cell distribution.This limits the applicability of the methodology developed in [8] and justifies the investigation presented in this work since it is a more general technique (with its corresponding advantages and disadvantages).Thus, we made a point of solving the calibration problem of the equal division rate by different techniques including Tikhonov-type regularization and Statistical methods.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 2015 estimation problem.Thus, in what follows, we shall consider the following simplified version of the model given by Equation (1):

Figure 1 :
Figure 1: Reconstructions of a non-smooth B using Tikhonov regularization (left), and statistical techniques (right).

Figure 2 :
Figure 2: Reconstructions of a smooth B using Tikhonov regularization (left), and statistical techniques (right).

Figure 4 :
Figure 4: The density N corresponding to the reconstructions of B using Tikhonov-type (Smoothing and Kullbacl-Leibler) regularization (left), and statistical techniques (right).