Mathematical modelling, analysis and numerical simulation for a general class of gene regulatory networks

Resumen

The research work developed in this thesis is mainly oriented to the mathematical modelling of biological systems, the behaviour of which is inherently stochastic, as it is the case of gene regulatory networks. Their relevance emerges from the fact that all necessary information for life cycle is encoded in the DNA. Consequently, the study of DNA expression, transcription into messenger RNA and translation into proteins, together with their regulation becomes essential to predict cells response to environmental signals. The inherent stochastic nature of gene expression makes these systems to be far away from the classical kinetic limit where the (macroscopic) deterministic methods are valid. In modelling these systems, we need to employ microscopic methods which incorporate the underlying stochastic behaviour. The Chemical Master Equation (CME) remains at the basis for the modelling of these phenomena. However, a closed form solution of the CME is unavailable in general, due to the large number (eventually infinity) of coupled equations. A widespread technique to approximate the CME solution is the Stochastic Simulation Algorithm (SSA), a computationally involved Monte Carlo type method. Although many numerical approximations emerge to reduce the complexity of the CME, we will focus on the Partial Integro-Differential Equation (PIDE) or Friedman model, which represents a continuous approximation of the CME. For the one dimensional version (self-regulation), the PIDE model has an analytic solution for its steady state. This fact will allow us to characterize the regions in the space of parameters in which the system changes its behaviour (unimodal, bimodal). Also we have carried out an stability analysis by means of entropy methods. Moreover, we obtain a multidimensional version of the Friedman model to handle more complex gene regulatory networks with more than one gene. The mathematical properties of the corresponding equation will be exhaustively analyzed, also with special emphasis on stability of the solution using entropy methods. In addition, we propose two semi-Lagrangian methods for the numerical solution of the multidimensional model. The first method results very efficient and scalable to higher dimensions, as the numerical results illustrate, although in practice exhibits first order convergence in time and space. Solutions provided by the proposed method are compared with those obtained by SSA to assess efficiency, accuracy and computational costs. For the second semi-Lagrangian method we develop the theoretical numerical analysis, thus proving second order convergence in time and space. This is clearly illustrated by a numerical example. However, the computational cost of this second approach results much higher, so that the scalability to higher dimensions seems a difficult task. All the numerical techniques have been implemented on a user friendly toolbox (SELANSI) which is detailed in the Appendix.