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

Abstract

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.

El trabajo de investigación se orienta al modelado matemático de sistemas biológicos que debido a su naturaleza presentan un carácter estocástico, entre los que se encuentran las redes de regulación genética. La importancia de estas radica en que toda la información para las funciones vitales de los seres vivos está codificada en el ADN. Por tanto, el estudio de la expresión del ADN, transcripción en ARN mensajero y traducción en proteinas, junto con su regulación, se convierte en esencial para predecir las respuestas de las células a se˜nales ambientales. La estocasticidad de los sistemas de regulación genética les hace estar lejos del clásico límite cinético para el cual los métodos deterministas (macroscópicos) tienen validez. En el modelado de estos sistemas surge la necesidad de emplear modelos microscópicos que incorporen el carácter estocástico subyacente. La Chemical Master Equation (CME) está en la base para el modelado de estos fenómenos. Sin embargo, no es posible obtener una solución de las CME en la mayoría de casos debido al gran número, incluso infinito, de ecuaciones acopladas. Una de las técnicas más extendidas para obtener esta solución es el empleo de métodos de Montecarlo, como el Stochastic Simulation Algoritm (SSA), aunque es numéricamente costoso. A pesar de la existencia de muchas aproximaciones para reducir la complejidad de la CME, nos centramos en las ecuaciones integro-diferenciales parciales o modelos de Friedman. Son una aproximación continua de las CME, que en su versión unidimensional (autorregulación) tienen una solución analítica para el estado estacionario. Este hecho nos permitirá caracterizar las regiones del espacio de parámetros en las cuales el sistema cambia su comportamiento (unimodal, bimodal). Además, hemos realizado un análisis de su estabilidad usando métodos de entropía. También se ha desarrollado una versión multidimensional del modelo de Friedman para poder manejar redes de regulación genética más complejas con más de un gen. Las propiedades matemáticas de este modelo serán analizadas en profundidad, haciendo especial hincapié en la estabilidad de la solución usando métodos de entropía. Además, se han propuesto dos métodos semi-Lagrangianos para la resolución numérica del modelo multidimensional desarrollado. El primer método resulta ser muy eficiente y escalable para dimensiones altas, como muestran los resultados numéricos, a pesar de tener orden uno de convergencia en espacio y tiempo. Las soluciones obtenidas con el método propuesto se comparan con aquellas obtenidas usando el SSA para evaluar su eficiencia, precisión y costes computacionales. Para el segundo método semi-Lagrangiano realizamos un análisis numérico teórico, probando orden de convergencia dos tanto en tiempo como en espacio. Esto se ilustrada claramente en un ejemplo numérico. Sin embargo, el coste computacional de este segundo método es mucho más alto, provocando que la escalabilidad para dimensiones altas sea una tarea difícil. Todas estas técnicas de resolución numérica han sido programadas e incluidas en una herramienta informática de fácil uso (SELANSI), que se detalla en el Apéndice.

O traballo de investigación oriéntase ó modelado matemático de sistemas biolóxicos que debido a súa natureza presentan un carácter estocástico, entre os que se atopan as redes de regulación xenética. A importancia destas radica en que toda a información para as funcións vitais dos seres vivos está codificada no ADN. Polo tanto, o estudo da expresión do ADN, transcrición en ARN mensaxeiro e tradución en proteinas, xunto coa súa regulación, convértese esencial para predicir as respostas das células a sinais ambientais. A estocasticidade dos sistemas de regulación xenética fainos estar lonxe do clásico límite cinético para o cal os métodos deterministas (macroscópicos) te˜nen validez. No modelado destes sistemas xurde a necesidade de empregar modelos microscópicos que incorporen o carácter estocástico subxacente. A Chemical Master Equation (CME) está na base para o modelado destes fenómenos. Sen embargo, non é posible obter unha solución das CME na maioría de casos debido ó gran número, incluso infinito, de ecuacións acopladas. Unha das técnicas máis estendidas para obter esta solución é o uso de métodos de Monte Carlo, como o Stochastic Simulation Algoritm (SSA) aínda que é numéricamente custoso. A pesar da existencia de moitas aproximacións para reducir a complexidade da CME, centrámonos nas ecuacións integro-diferenciais parciais ou modelos de Friedman. Son unha aproximación continua das CME que na súa versión unidimensional (autorregulación) te˜nen unha solución analítica para o estado estacionario. Este feito permitiranos caracterizar as rexións no espazo de parámetros nas que o sistema cambia o seu comportamento (unimodal, bimodal). Ademais levouse a cabo unha análise da súa estabilidade usando métodos de entropía. Desenvolveuse unha versión multidimensional do modelo de Friedman para poder manexar redes de regulación xenética máis complexas con máis dun xene. As propiedades matemáticas deste modelo serán analizadas en profundidade, facendo especial fincapé na estabilidade da solución usando métodos de entropía. Ademais, propuxéronse dous métodos semi-Lagranxianos para a resolución numérica do modelo multidimensional desenvolto. O primeiro método resulta ser moi eficiente e escalable para dimensións altas, como mostran os resultados numéricos, a pesares de ter orde un de converxencia en espazo e tempo. As solucións obtidas co método proposto compáranse con aquelas obtidas usando o SSA para avaliar a súa eficiencia, precisión e custes de computación. Para o segundo método semiLagranxiano realizamos unha análise numérica teórica, probando orde dúas de converxencia en tempo e espazo. Isto é claramente ilustrado nun exemplo numérico. Sen embargo, o custe de computación deste segundo método é moito máis alto provocando que a escalabilidade para dimensións altas sexa unha tarefa difícil. Todas estas técnicas de resolución numérica foron programadas e incluidas nunha ferramenta informática de fácil uso (SELANSI), que é detallada no Apéndice.