German Title: Numerische Methoden zur Parameterschätzung in Dynamischen Systemen mit Rauschen mit Anwendungen in der Systembiologie

Preview

PDF, English - main document Download (20MB) | Terms of use

Abstract

This thesis comprises the modelling of and parameter estimation in dynamical systems, with a focus on applications in systems biology. In an interdisciplinary research project on the systems biology of cancer, we develop a predictive mathematical model of an intracellular crosstalk in cytokine signalling. Expected and unexpected predictions are confirmed in experiments and lead to new biological insights. For model calibration with measurement data, we apply well established methods for parameter estimation in ordinary differential equation models. Extending these to stochastic differential equations, we develop, analyse, and implement a new method for parameter estimation in dynamical processes with noise, and demonstrate its performance in several selected examples from systems biology and mathematical finance.

Many processes, especially in biology, obey deterministic ground rules (e.g. metabolic processes or signal transduction pathways), but may be heavily influenced by fluctuations and stochasticity inherent to the system that change its behaviour both qualitatively and quantitatively. Therefore, frequently, a deterministic description is not constructive.

A large class of such systems can be adequately described by nonlinear multi-dimensional stochastic differential equations (SDEs). Classical estimation techniques for SDEs, relying on (approximations of) transition densities, are all too often not applicable to these problems, due, inter alia, to their high computational costs and prerequisites on the measurements. The proposed new method is based on the method of multiple shooting, using piecewise deterministic solutions of ordinary differential equations (ODEs) to approximate the SDE realization that corresponds to the studied process from which measurements have been taken.

The generally discontinuous concatenation of ODE trajectories mimics the consequences of stochastic effects, and, further, allows to formulate the parameter estimation problem as a deterministic nonlinear optimization problem that can be solved with efficient derivative-based solution methods. In this thesis, a generalized GAUSS-NEWTON method is deployed.

Main results and contributions of this thesis are summarized in the following:

• We propose a new method for parameter estimation in nonlinear multi-dimensional SDEs, based on a piecewise approximation by solutions of ODEs. Discontinuities (jumps) occurring at the interval borders are used for regularization. Unknown parameters and initial states are estimated by a generalized weighted least squares method from data that can originate from direct complete or partial state measurements or from indirectly observed quantities. Measurement data may be afflicted with errors and arbitrarily sampled. Non-linear parameter and point constraints may be formulated as equality and inequality constraints. The resulting nonlinear constrained optimization problems are highly structured and efficiently solved using a generalized GAUSS-NEWTON method.

• We give a proof that the discontinuities at the interval borders asymptotically tend to zero if the number of equidistantly distributed shooting nodes goes to infinity.

• We show in a numerical analysis that the resulting equation systems are sparse, that the number of nonzero elements depends only linearly on the number of shooting nodes, and give sharp upper bounds. Moreover, we prove that the sparsity is maintained if an appropriate (stable) decomposition is a applied.

• It is demonstrated in comparative simulation studies that the estimates are robust w.r.t. the exact choice of jump regularization weights. Moreover, the effects of jump regularization on estimates and approximated trajectories are investigated and described.

• A lifting approach with per-interval parameter sets, coupled by additional equality constraints, is developed and its numerical properties are analysed. Moreover, we propose a homotopy method for the treatment of hard problems.

• We demonstrate the performance of the new estimation technique in examples from systems biology, each shedding some light on different aspects. Especially, we show that the method can also be used for hidden state estimation and for trajectory reconstruction in time spans without observations. Further, we derive a criterion for local grid refinement.

• We show for an ORNSTEIN-UHLENBECK process driven by a LEVY jump process, that, in addition to mean reversion level and mean reversion rate, also the diffusion constant may be estimated by analysing the jump residuals.

• The software package :sfit is an efficient implementation of the proposed method, offering easy symbolic problem formulation to the user, from which the stochastic parameter estimation problem can be automatically built and solved.

Parts of this work emerged in the interdisciplinary research project SBCancer of the Helmholtz Alliance on Systems Biology. In close collaboration with expert biologists, we developed an mathematical model for a crosstalk of two cytokines in human skin cells that interfere in a signalling pathway frequently found aberrantly activated in cancer. After an extensive analysis of the deployed measurement data processing, the model proposed by the author of this thesis has been calibrated from experimental data. Its counter-intuitive predictions have been verified in wet lab experiments and lead to new biological insights.

Main novelties and contributions in this thesis are:

• Development of a mathematical model of the crosstalk of two cytokines in human keratinocytes (HaCaT cell line). The predicted and hitherto unknown nonlinear moderating effects of GM-CSF on the IL-6-induced JAK-STAT signalling pathway has been verified in vitro.

• An extensive mathematical analysis of the frequently utilized quantitative WESTERN blotting measurement procedure shows that established data normalization methods, relying on housekeeping proteins or manually added calibrator proteins, are prone to signal deteriorating statistical artefacts. Moreover, we show that the frequently declared assumption of normally distributed measurement errors cannot be maintained if these normalization techniques are applied.

• We propose as a remedy a normalization technique based on the calculation of amplification factors, and develop criteria for (approximate) normally distributed errors. These criteria can be easily checked using solely the raw measurement data. Moreover, we demonstrate the advantages of the proposed amplification factors method in a large comparative simulation study.

Translation of abstract (German)

Die vorliegende Arbeit behandelt die Modellierung von und Parameterschätzung in dynamischen Systemen mit Fokus auf Anwendungen in der Systembiologie. In dieser interdisziplinären Forschungsarbeit zur Systembiologie von Krebserkrankungen wird ein prädiktives mathematisches Modell intrazellulärer Wechselwirkungen zweier Zytokinsignaltransduktionswege entwickelt und im Experiment bestätigt. Unerwartete Vorhersagen führen auf neue biologische Erkenntnisse. Zur Modellkalibrierung kommen klassische Methoden der Parameterschätzung in gewöhnlichen Differentialgleichungen zum Einsatz. Daran anknüpfend wird in dieser Arbeit eine neue numerische Methode zur Parameterschätzung in stochastischen Einflüssen unterliegenden dynamischen Prozessen entwickelt, theoretisch untersucht, implementiert und ihre Leistungsfähigkeit an Beispielen aus Systembiologie und Finanzwissenschaft demonstriert.

Viele Prozesse, insbesondere biologische, zeigen ein grundsätzlich gerichtetes Verhalten, das zwar gewissen Gesetzmäßigkeiten unterliegt (etwa Stoffwechselprozesse oder Signaltransduktionen), allerdings auch durch intrinsische Zufälligkeiten qualitativ und quantitativ erheblich beeinflusst wird, sodass eine rein deterministisch-mechanistische Modellierung der auftreten den Vorgänge oft nicht zielführend ist.

Eine große Klasse solcher Systeme lässt sich durch nichtlineare mehrdimensionale stochastische Differentialgleichungen (SDEs) adäquat beschreiben. Klassische stochastische Schätzverfahren, die auf (Approximation von) ¨Ubergangswahrscheinlichkeiten basieren, sind auf diese Problemklasse u. a. wegen ihres Rechenaufwands und speziellen Anforderungen an die Messdaten oft nicht oder nur eingeschränkt anwendbar.

In dieser Arbeit wird ein neues Mehrzielverfahren entwickelt, das die dem beobachteten Prozess entsprechende Realisierung der SDE stückweise durch deterministische Lösungen gewöhnlicher Differentialgleichungen (ODEs) approximiert. Diese i. d. R. unstetige Konkatenation der Trajektorien erlaubt zum einen die Wiedergabe stochastischer Effekte, zum anderen die Darstellung des Parameterschätzproblems als deterministisches nichtlineares Optimierungsproblem, das unter Verwendung effizienter ableitungsbasierter Verfahren gelöst werden kann. In der vorliegenden Arbeit kommt hierzu ein verallgemeinertes GAUSS-NEWTON-Verfahren zum Einsatz.

Die wesentlichen Ergebnisse und Resultate der vorliegenden Arbeit umfassen insbesondere:

• Es wird ein neues Verfahren zur Parameterschätzung in mehrdimensionalen nichtlinearen SDEs entwickelt, das auf stückweiser Approximation durch Lösungen von ODEs basiert. An Intervallgrenzen auftretende Sprünge werden zur Regularisierung verwendet. Unbekannte Parameter und initiale Systemzustände werden mittels einer verallgemeinerten, gewichteten Kleinste-Quadrate-Methode aus Messdaten geschätzt, die aus direkten Zustandsmessungen oder aus durch Messfunktionen beschriebenen indirekten Beobachtungen stammen, mit Fehlern behaftet und an beliebigen Zeitpunkten erfasst worden sein können. Nichtlineare Parameter- und punktweise Zustandsbeschränkungen können als Gleichungs- und Ungleichungsbedingungen formuliert werden. Die resultierenden nichtlinearen beschränkten Optimierungsprobleme werden unter Strukturausnutzung mit einem verallgemeinerten GAUSS-NEWTON-Verfahren effizient gelöst.

• Es wird ein Beweis gegeben, dass im Falle einer asymptotisch gegen unendlich gehenden Zahl von äquidistanten Mehrzielknoten die Unstetigkeiten (Sprünge) an den Intervallgrenzen gegen null gehen.

• Eine numerische Analyse offenbart die Dünnbesetztheitsstruktur der auftretenden Gleichungssysteme. Es wird gezeigt, dass die Anzahl der Nichtnulleinträge nur linear in der Zahl der Mehrzielknoten (Zeitgitter) wächst und einer scharfen oberen Schranke genügt. Zusätzlich wird bewiesen, dass sich unter Verwendung einer geeigneten (stabilen) Zerlegung die Dünnbesetztheit erhalten lässt.

• Eine vergleichende Simulationsstudie zeigt, dass die Schätzungen sehr robust hinsichtlich der Wahl der Sprunggewichte sind. Zusätzlich wird dargelegt, welchen Einfluss die Intensität der Sprungregularisierung auf Schätzungen und approximierte Trajektorien hat.

• Es wird ein Lifting-Ansatz mit intervallweisen Parametersätzen, gekoppelt durch zusätzliche Gleichungsbedingungen, entwickelt und seine numerischen Eigenschaften analysiert. Weiter wird ein Homotopieverfahren zur Behandlung schwer lösbarer Probleme vorgeschlagen.

• Die Leistungsfähigkeit der Methode wird an Beispielen aus der Systembiologie, die jeweils unterschiedliche Aspekte beleuchten, demonstriert. Insbesondere wird gezeigt, dass sich die Methode auch zur Rekonstruktion von unbeobachteten Zuständen (engl. hidden state estimation), sowie zur vollständigen Rekonstruktion von Trajektorien in Zeiträumen, in denen keine Beobachtungen vorliegen, einsetzen lässt. Weiter wird ein Kriterium zur lokalen Verfeinerung des Mehrzielgitters gegeben.

• Für einen durch einen LEVY-Sprungprozess getriebenen ORNSTEIN-UHLENBECK-Prozess (modifiziertes VASICEK-Zinsmodell) wird gezeigt, dass sich neben Gleichgewichtsniveau und Steifigkeit aus den Sprungresiduen auch die Diffusionskonstante schätzen lässt.

• Im Softwarepaket :sfit wird eine effiziente Implementierung der in dieser Arbeit entwickelten Parameterschätzmethode bereitgestellt, die auf Nutzerseite eine einfache symbolische Problemformulierung erlaubt.

Teile dieser Arbeit entstanden im interdisziplinären Forschungsprojekt SBCancer der Helmholtz-Allianz Systembiologie und beschäftigen sich mit der Modellierung der Wechselwirkung zweier an Krebs beteiligter Zytokin-Signalkaskaden in menschlichen Hautzellen. Vom Autor wurde in enger Kooperation mit den beteiligten Biologen ein mathematisches Modell entwickelt und durch Messdaten kalibriert. Die mittels diesem Modell berechneten kontraintuitiven Vorhersagen eines vorgeschlagenen Doppelstimulationsexperimentes konnten im Labor bestätigt werden und führten auf neue biologische Erkenntnisse.

Wesentliche Beiträge und Neuerungen in dieser Arbeit sind die folgenden:

• Entwicklung eines mathematischen Modells eines Crosstalks der Signalwege zweier Zytokine in menschlichen Keratinozyten (Zelllinie HaCaT). Die vom Modell vorhergesagte und zuvor unbekannte nichtlineare moderierende Wirkung von GM-CSF auf den IL-6-induzierten JAK-STAT-Signalweg wurde im Labor in vitro nachgewiesen.

• Eine aufwendige mathematische Analyse des in der Zellbiologie häufig eingesetzten quantitativen WESTERN-blot-Messverfahrens zeigt, dass etablierte Normalisierungstechniken, die auf in angenommen konstanter Konzentration vorliegenden Haushaltsproteinen (engl. house-keeping proteins) oder auf manuell hinzugefügten Kalibrierungsproteinen beruhen, extrem anfällig für signalzerstörende statistische Artefakte sind. Zusätzlich wird gezeigt, dass sich die häufig getroffene Annahme normalverteilter Fehler nach Anwendung dieser Normalisierung i. d. R. nicht halten lässt.

• Eine auf Verstärkungsfaktoren beruhende Methode wird als Ersatz vorgeschlagen. Es werden leicht nachzuprüfende Kriterien an die Rohdaten entwickelt, anhand derer entschieden werden kann, ob die auf diese Weise prozessierten Messdaten weiterhin (approximativ) normalverteilt sind. Des Weiteren werden in einer vergleichenden Simulationsstudie die Vorteile der vorgeschlagenen Verstärkungsfaktormethode (engl. amplifications factors method) demonstriert.