Section: New Results

Estimation for complex and biological systems

Participants: T. Bastogne, C. Lacaux, S. Mézières, S. Tindel, P. Vallois

Tumor growth modeling

This project is an extension of our article [15] , which will be written in 2014. A cancer tumor can be represented for simplicity as an aggregate of cancer cells, each cell behaving according to the same discrete model and independently of the others. Therefore to measure its size evolution, it seems natural to use tools coming from dynamics of population, for instance the logistic model. This deterministic framework is well-known but the stochastic one is worthy of interest. We are currently working on a model in which we suppose that the size $V_t$ at time $t$ of the tumor is a diffusion process of the type :

where ${\displaystyle {\left(B_t\right)}_{t\ge 0}}$ is a standard brownian motion starting from zero.Then (i) We define a family of time continuous Markov chains which models the evolution of the rate of malignant cells and approximate (under some conditions) the diffusion process $\left(V_t\right)$. (ii) We study in depth the solution to equation (1 ). This diffusion process lives between two frontiers : 0 and $\kappa$. In this stochastic setting, the role of $\kappa$ is not so clear and we contribute to understand it. We describe the asymptotic behavior of the diffusion according to the values of the parameters. The tools we resort to are boundary classification criteria and Laplace transform of the hitting time to biological worthwhile level. We believe we are able in particular to express the mean of the hitting time.

The next step in this project can be summarized as follows: at this point in our investigations on tumor growth modelization, we have identified a pertinent and consistent model. Nevertheless our study remains theoretical. A statistical estimation of the parameters $r,\kappa ,c$ $\beta$ is thus in order. This would permit to apply our model to real data. A further objective could be to consider a more complex form for the logistic term, see e.g. Schurtz (2007).

Local score associated with long biological sequences

Statistical properties of the distribution of the local score is largely used by molecular biologists to extract important features in biological sequences and in particular to determine the most significant one among a collection of biological sequences. The probabilistic model which is commonly used is the following. Associated with a sequence ${\left({ϵ}_i\right)}_{i\ge 1}$ of independent, centered and reduced random variables, consider $S_n={ϵ}_1+\cdots +{ϵ}_n$ and

In biological sequence analysis, $\left({ϵ}_i\right)$ can for example correspond to the physical or chemical properties of the $i$-th amino acid or nucleotid of the sequence ; it can also reflect the similarity between components of two sequences. The local score ${\overline{U}}_n$ is the supremum of $\left(U_n\right)$ up to time $n$. Molecular biologists are interested by this supremum as it highlights the best part of the studied sequence, the eventual segment of DNA transmitted by a common ancestor for sequence comparison or the best hydrophobic segment of a protein that would thus naturally move in a transmembrane place. It is clear that the trajectory of $\left(U_n\right)$ can be decomposed in a succession of 0 and excursions above 0. These excursions have an important biological interpretation and in particular the highest one corresponds to the best segment due to the physico chemical property or similarity scores that have been chosen. Note that the local score ${\overline{U}}_n$ can be viewed as the maximum of the heights of all the excursions up to time $n$. In the article [22] , we are interested in complete excursions up to a fixed time. This leads us to introduce the maximum $U_n^{*}$ of the heights of all the excursions up to time $n$. The second variable which will play an important role is ${\theta }_n^{*}$ the time necessary to reach its maximal height $U_n^{*}$.

We believe that the knowledge of the joint distribution of the pair $(U_n^{*},{\theta }_n^{*})$ would permit to get more efficient statistical tests than the ones only based on the local score. This point should be developed in a forthcoming paper.

However, it seems difficult to determine explicitly the law of $(U_n^{*},{\theta }_n^{*})$ for a fixed $n$. This difficulty can be overcome considering biological sequences which have a large number of bases and approximating the initial random walk $\left(S_n\right)$ by a Brownian motion $\left(B_t\right)$ started at 0. Using the functional theorem of convergence of Donsker, the process $\left(U_k\right)$ can be compared to

This leads us to consider:

1. the local score $\overline{U}\left(t\right)$ which is the maximum of the heights of all the excursions of $U\left(s\right)$ up to time $t$,

2. the maximum $U^{*}\left(t\right)$ of the heights of all the complete excursions up to time $t$,

3. the time ${\theta }^{*}\left(t\right)$ taken by $U\left(s\right)$ starting from the beginning of the largest excursion to hit the maximal level $U^{*}\left(t\right)$.

The approximation of $\left(U_n\right)$ by $\left({\widehat{U}}_t\right)$ permits to prove that the asymptotic distribution of ${\displaystyle \left(\frac{U_n^{*}}{\sqrt{n}},\frac{{\theta }_n^{*}}{n}\right)}$ as $n\to \infty$ is the one of $\left(U^{*}\left(1\right),{\theta }^{*}\left(1\right)\right)$. Consequently, our initial problem in the discrete setting reduces to determine the joint law of $(U^{*}\left(t\right),{\theta }^{*}\left(t\right))$, where $t>0$ is given. We determine in [22] the distribution and the density functions of $\left(U^{*}\left(t\right),{\theta }^{*}\left(t\right)\right)$.

Bacteriophage therapies

In the last years Bacteriophage therapies are attracting the attention of several scientific studies. They can be a new and powerful tool to treat bacterial infections or to prevent them applying the treatment to animals such as poultry or swine. Very roughly speaking, they consist in inoculating a (benign) virus in order to kill the bacteria known to be responsible of a certain disease. This kind of treatment is known since the beginning of the 20th century, but has been in disuse in the Western world, erased by antibiotic therapies. However, a small activity in this domain has survived in the USSR, and it is now re-emerging (at least at an experimental level). Among the reasons of this re-emersion we can find the progressive slowdown in antibiotic efficiency (antibiotic resistance). Reported recent experiments include animal diseases like hemorrhagic septicemia in cattle or atrophic rhinitis in swine, and a need for suitable mathematical models is now expressed by the community.

At a mathematical level, whenever the mobility of the different biological actors is high enough, bacteriophage systems can be modeled by a kind of predator-prey equation. Namely, set $S_t$ (resp. $Q_t$) for the bacteria (resp. bacteriophages) concentration at time $t$. Then a model for the evolution of the couple $(S,Q)$ is as follows:

where $\alpha$ is the reproducing rate of the bacteria and $k$ is the adsorption rate. In equation (3 ), $d$ also stands for the quantity of bacteriophages inoculated per unit of time, $m$ is their death rate, we denote by $b$ the number of bacteriophages which is released after replication within the bacteria cell, $\zeta$ is the delay necessary to the reproduction of bacteriophages (called latency time) and the coefficient $e^{-\mu \zeta }$ represents an attenuation in the release of bacteriophages (given by the expected number of bacteria cell's deaths during the latency time, where $\mu$ is the bacteria's death rate). A given initial condition $(S_0,Q_0)$ is also specified, and the term $\epsilon \mathrm{d}{W}_t$ takes into account a small external noise standing for both uncertainties on the measures and the experiment conditions. One should be aware of the fact that the latency time $\zeta$ (which can be seen as the reproduction time of the bacteriophages within the bacteria) cannot be neglected, and is generally of the same order (about 20mn) as the experiment length (about 60mn).

With this model in hand, our main results in this direction (see [1] ) have been the following:

• Quantification of the exponential convergence to a bacteria-free equilibrium of equation (3 ) when $d$ is large enough.

• Use of the previous result plus concentration inequalities in order to study the convergence of the noisy system to equilibrium in a reasonable time range.

• Simulation of the stochastic processes at stake in order to observe the convergence to equilibrium.

Light transport in tissues with probabilistic methods

Photodynamic therapy (PDT) is a type of phototherapy used for treating several diseases such as acne, bacterial infection, viruses and some cancers. The aim of this treatment is to kill pathological cells with a photosensitive drug that is absorbed by the target cells and that is then activated by light. For appropriate wavelength and power, the light beam makes the photosensitizer produce singlet oxygen at high doses and induces the apoptosis and necrosis of the malignant cells. Our project focuses on an innovative application: the interstitial PDT for the treatment of high-grade brain tumors. This strategy requires the installation of optical fibers to deliver the light directly into the tumor tissue to be treated, while nanoparticles are used to carry the photosensitizer into the cancer cells. In order to optimize the intra-cerebral position of our optical fiber, two fundamental questions have to be answered:

1. What is the optimal shape and position of the light source in order to optimize the damage on malignant cells?

2. Is there a way to identify the physical parameters of the tissue which drive the light propagation?

Notice that we are obviously not the first ones to address these issues, and there is nowadays a consensus in favor of the algorithm proposed by L. Wang and S. L. Jacques for the simulation of light transport in biological tissues. However, our starting point is the observation that the usual methods slightly lack of formalism and miss formal representations that answer the questions of identifiability. In [25] , in the framework of homogeneous biological tissues, we propose an alternative MC method to Wang's algorithm. Then we also propose a variance reduction method. Interestingly enough, our formulation also allows us to design quite easily a Markov chain Monte Carlo (MCMC) method based on Metropolis-Hastings algorithm and to handle the inverse problem (of crucial importance for practitioners), consisting in estimating the optical coefficients of the tissue according to a series of measurements. We have compared the proposed MC and MCMC method and Wang's algorithm: we see that our MC method is much more consistent. However, MCMC methods induce quick mutations, which paves the way to very promising algorithms in the inhomogenous case. To handle the inverse problem, we derive a probabilistic representation of the variation of the fluence with respect to the absorption and scattering coefficients. This leads us to the implementation of a Levenberg-Marquardt type algorithm that gives an approximate solution to the inverse problem.

System identification of gap junctional intercellular communication channels of two cancer cell lines.

The main challenge addressed in this work [12] , [14] was to assess the relevance of a proposed model-based approach to detect differences between gap junctional intercellular communication channels of two cancer cell lines. To that aim, KB and FaDu, two human head and neck carcinoma cell lines, were used. The former expresses connexin proteins (positive line) while the latter does not (negative line). Moreover, each cell line was tested on spheroid (3-D) and monolayer (2-D) slices and in vitro assays were repeated six times. Continuous-time system identification algorithms of the Matlab System Identification and CONTSID toolboxes are tested and applied to a set of in vitro data. Results firstly show an acceptable fit of the biological responses but they mainly emphasize the possibility to use several model parameters as statistics to discriminate different cancer cell lines. So, this study exemplifies the potential contribution of dynamic system identification methods and tools to the discovery of new diagnostic biomarkers in systems biology.

Photodynamic therapy modeling and control.

We have also carried on the development of methodological and technological innovations for the realtime control of the therapeutic efficiency in photodynamic therapy (Tylcz:2013). One part of the innovation has been protected by a patent submitted in 2012 (No.1261339, INPI) and extended in 2013. A demonstration platform is currently in development.

Bio-inspired system reliability method.

Based on previously developed works (Keinj, 2011, 2012), we have also proposed in [15] an extension of the target theory in biology applied to system reliability. In this development, we consider rough products produced by a factory. Each product coming from the plant has $m$ vital elements and some elements can be damaged. To obtain a perfect product (i.e. all the constitutive $m$ elements are safe) all the damaged elements are repaired and a test phase follows. The result of this two-steps procedure is random. We suppose that the number $Z_k$ of non-damaged elements is a Markov chain valued in the set $\{0,1,\cdots ,m\}$, where $k$ is the number of applied repairing-test phases. We have a qualitative result which says that if the repair phase is efficient then $P(Z_k=m)$ is close to 1. As for production of a large number $n$ of products, the former result allows us to give conditions under which either the $n$ elements or a fraction of these $n$ elements are (is) safe after the application of $k$ previous maintenance phases.

Dynamical Global Sensitivity Analysis as an Early Warning for System's Critical Transition.

In biology, systems with bifurcations may experience abrupt irreversible and often unwanted shifts in their performance, called critical transitions. For many systems like climate, economy, ecosystems it is highly desirable to identify indicators serving as early warnings of such regime shifts. Several statistical measures were recently proposed as early warnings of critical transitions including increased variance, autocorrelation and skewness of experimental or model-generated data. The lack of automatized tool for model-based prediction of critical transitions led to designing DyGloSA, a Matlab toolbox for dynamical global parameter sensitivity analysis (GPSA) of ordinary differential equations models. One part of our research activity in 2013 was focused on the implementation of a global sensitivity analysis method developed in (Dobre, 2011, 2012) into DyGloSA for dynamical global parameter sensitivity analysis (GPSA) of ordinary differential equations models. This work has been carried out in the context of a collaboration with the University of Luxembourg and more precisely the Thomas Sauter's team. We have shown in [2] that tools developed in this toolbox are efficient to analyze several models with bifurcations and predict the time periods when systems can still avoid going to a critical transition by manipulating certain parameter values, which is not detectable with the existing SA techniques.