Parameteridentifikation
(Bioprozess-Modell)

(Stand: Juni 2022)

Die vorgestellte Anwendung demonstriert am Beispiel eines simplifizierten UpStream-Prozesses zur Herstellung eines rekombinanten Bioproduktes eine Strategie zur Parameter-Identifikation. Ziel der Parameter-Identifikation ist die Parametrierung eines mathematischen Bioprozess-Modells, welches das essentielle System-Verhalten des Bioprozesses abbildet. Im vorliegenden Fall erfolgt die Bildung des Bioproduktes mithilfe eines genetisch veränderten Organismus. Das parametrierte Bioprozess-Modell kann im nachfolgenden zur Analyse (z. B. Sensitivitätsanalyse) oder zur Optimierung (z. B. Ausbeute-Maximierung) des bestehenden Bioprozesses eingesetzt werden.

Anwendungsbeispiel


Die Kultivierung des genetisch veränderten Organismus $\color{#4daf4a}{X}$ (vgl. E. coli) erfolgt in einem ideal durchmischten klassischen Rührkessel-Reaktor. Die Bildung des rekombinanten Bioproduktes $\color{#377eb8}{P}$ wird durch die Zugabe eine Aktivators $a$ (vgl. IPTG) zum Zeit­punkt $t_{\text{a}}$ induziert [1]. Hierbei konkurrieren Bio­mas­se­wachs­tum und Pro­dukt­bil­dung je­weils um dasselbe Substrat $\color{#ff7f00}{S}$ (vgl. Glu­ko­se). Bei Ablauf der Pro­zess­lauf­zeit $t_{\text{harvest}}$ erfolgt die Ernte und Aufarbeitung des gebildeten Bioproduktes. Desweiteren sei bekannt, dass hohe Substrat-Konzentrationen $\color{#ff7f00}{S}$ das Biomassewachstum hemmen. Nach Zugabe des Aktivators wird ein Großteil des Substrates durch Produktbildung verbraucht. Zusätzlich sei eine deutliche Inhibition des Biomassewachstums durch die ansteigende Produktkonzentration gegeben. Nach Auszehrung des Substrates stagniert die Produktbildung. Ebenfalls kann der Zerfall von Biomasse beobachtet werden.

Voraussetzungen zur Durchführung der Parameter-Identifikation sind eine ausreichende experimentelle Datenbasis sowie ein entsprechendes mathematisches Bioprozess-Modell. Die Datenbasis muss das betrachtete System mit mindestens allen essentiellen Eigenschaften abbilden. Hierzu sind eine Vielzahl von Experimenten bzw. Versuchen durchzuführen, um Messdaten für die Parameter-Identifikation zu gewinnen. In den Experimenten werden festgelegte Steuergrößen in geeigneter Kombination variiert. Steuergrößen sind jene Größen, welche einen Eingriff in das jeweilige System erlauben. Im betrachteten Beispiel werden die Start-Biomassekonzentration $X_0$, die Start-Substratkonzentration $S_0$, der Aktivierungszeitpunkt der Produktbildung $t_{\text{a}}$ sowie der Ernte-Zeitpunkt $t_{\text{harvest}}$ als Steuergrößen festgelegt.

Mittels Versuchsplanung inklusive der Aussteuerung von Steuergrößen können relevante Experimente geplant und durchgeführt werden. Im betrachteten Anwendungsbeispiel sei die Datenbasis durch 9 Experimente gegeben, welche in nachfolgender Abbildung visualisiert sind (zum Erweitern anklicken):

Messdaten

Desweiteren sei das notwendige mathematische Bioprozess-Modell bekannt und in Form des nachfolgenden Differentialgleichungssystems gegeben [2,3]:

\begin{align} \dot{X}=&\,\mu(S,P)\cdot X -k_{\text{d}}\cdot X \label{eq:ode_cX} \\[10pt] \dot{S}=&\,-\dfrac{\mu(S,P)\cdot X}{Y_{\text{X,S}}} -\dfrac{\rho(S,a)\cdot X}{Y_{\text{P,S}}} - m \cdot X \label{eq:ode_cS} \\[10pt] \dot{P}=&\,\rho(S,a)\cdot X \label{eq:ode_cP} \end{align}

Die Parameter $k_{\text{d}}$ und $m$ entsprechen hierbei den konstanten Raten des Zerfalls von Biomasse sowie dem Verbrauch von Substrat zur Bedienung des Erhaltungsstoffwechsels der Zelle. Die beiden Ertragskoeffizienten $Y_{\text{X,S}}$ und $Y_{\text{P,S}}$ spiegeln das Verhältnis zwischen verbrauchtem Substrat zu hieraus gebildeter Biomasse bzw. gebildetem Produkt wider. Die nachfolgenden Gleichungen \eqref{eq:mu_kinetics} und Gl. \eqref{eq:pi_kinetics} bilden die Modell-Kinetiken für die Beschreibung von Biomasse-Wachstum bzw. -Zerfall und Produktbildung ab [2]:

\begin{eqnarray} \mu(S,P) & = & \frac{\mu_{\max}\cdot S}{K_{\text{S}}+S+\frac{S^2}{K_{\text{I,S}}}+\frac{S\cdot P}{K_{\text{I,P}}}}\label{eq:mu_kinetics}\\[10pt] \rho(S,a) & = & \begin{cases} 0 & \text{for }a=0\\ \frac{\rho_{\max}\cdot S}{K_{\text{P}}+S} & \text{for }a=1 \end{cases}\label{eq:pi_kinetics} \end{eqnarray}

Methodischer Hintergrund


Die Parameter-Identifikation, also die Bestimmung der Beträge aller Modell-Parameter, wird als mathematische Minimierung implementiert. Ziel ist es, diese Beträge so zu bestimmen, dass die Abweichung (Fehler) zwischen experimentellen Daten und simulierten Verläufen minimal wird. Hierzu stehen eine Vielzahl von etablierten Verfahren bereit. Im vorliegenden Anwendungsbeispiel wird die Methode der kleinsten Fehlerquadrate (engl.: least squares method) zusammen mit einem direkten heuristischen Suchverfahren (downhill-simplex [4]) eingesetzt. Desweiteren sei gegeben, dass alle Modell-Parameter strukturell sowie praktisch Identifizierbar sind [5].

Das folgende Schema visualisiert den Ablauf der Parameter-Identifikation und ihre wesentlichen Komponenten, deren Zusammenspiel danach beschrieben wird:

Ablauf Parameteridentifikation

Zur Durchführung der Parameter-Identifikation müssen Wertebereiche und Start-Werte $\theta_0$ für die (vektorisierten) Modell-Parameter $\theta$ vorgegeben werden. In den meisten Anwendungsfällen können die Begrenzungen (engl.: constraints oder bounds) aufgrund von Vorwissen sowie formal-kinetischen bzw. mathematischen Gegebenheiten eingeschränkt werden. Nachfolgende Tabelle listet die in diesem Anwendungsbeispiel verwendeten Modell-Parameter sowie deren entsprechende Wertebereiche (Begrenzungen) auf:

Symbol Bezeichnung Begrenzungen Einheit
untere obere
$\mu_{\max}$ maximale spezifische Wachstumsrate $0.1$ $0.7$ $\text{h}^{-1}$
$\rho_{\max}$ maximale spezifische Produktbildungssrate $0.01$ $0.05$ $\text{h}^{-1}$
$K_{\text{S}}$ Halbsättigungskonstante Substrat $0.5$ $5$ $\text{g}_{\text{S}} \cdot \text{L}^{-1}$
$K_{\text{P}}$ Halbsättigungskonstante Produkt $0.5$ $5$ $\text{g}_{\text{S}} \cdot \text{L}^{-1}$
$K_{\text{I,S}}$ Inhibitionskonstante Substrat $5$ $50$ $\text{g}_{\text{S}} \cdot \text{L}^{-1}$
$K_{\text{I,P}}$ Inhibitionskonstante Produkt $0.05$ $2$ $\text{g}_{\text{P}} \cdot \text{L}^{-1}$
$k_{\text{d}}$ spezifische Biomasse-Zerfallsrate $0$ $0.1$ $\text{h}^{-1}$
$m$ spezifische Substrat-Maintenance-Rate $0$ $0.1$ $\text{g}_{\text{S}} \cdot \text{g}_{\text{X}}^{-1} \cdot \text{h}^{-1}$
$Y_{\text{X,S}}$ Ertragskoeffizient Biomasse/Substrat $0.1$ $1$ $\text{g}_{\text{X}} \cdot \text{g}_{\text{S}}^{-1}$
$Y_{\text{P,S}}$ Ertragskoeffizient Produkt/Substrat $0.01$ $0.1$ $\text{g}_{\text{P}} \cdot \text{g}_{\text{S}}^{-1}$

Die Festlegung geeigneter Start-Werte stellt im allgemeinen eine große Herausforderung dar. Diese Start-Werte sind im Idealfall durch bestehendes Vorwissen bereits gegeben, oder können bei einfachen Modellen, mit sehr wenigen Parametern, empirisch ermittelt werden. Handelt es sich, wie im vorliegenden Anwendungsbeispiel, um ein dynamisches nicht-lineares Modell mit einer Vielzahl von Modell-Parametern, werden geeignete Start-Werte $\theta_0$ mithilfe verschiedener Methoden ermittelt. Die Generierung dieser Start-Wert-Kombinationen kann beispielsweise (voll-)kombinatorisch oder mittels verschiedener stochastischen bzw. probabilistischen Sampling-Strategien [6–8] erfolgen.

Die optimalen Beträge der Modell-Parameter $\Omega$ werden mittels Minimierung (Gl. \ref{eq:optimisation}) nachfolgender Zielfunktion $J\left(\theta\right)$ (Gl. \ref{eq:costfun}) ermittelt.

\begin{equation} \Omega=\underset{\theta}{\mathrm{argmin}}\left(J\left(\theta\right)\right) \label{eq:optimisation} \end{equation} \begin{equation} J\left(\theta\right)=\sum_{j=1}^{k}\sum_{i=1}^{n}\left(y_\text{data}\left(t_{j,i}\right)-y_\text{model}\left(\theta,t_{j,i}\right)\right)^{2} \label{eq:costfun} \end{equation}

Diese Zielfunktion summiert die quadrierten Differenzen zwischen der beobachteten Größe $y_{\text{data}}\left(t_{j,i}\right)$ und der simulierten Größe $y_{\text{model}}\left(\theta,t_{j,i}\right)$. Die Indices $i$ ($i=1,\ldots,n$) und $j$ ($j=1,\ldots,k$) referenzieren den diskreten Wert ($i$) innerhalb der Zeitreihe des jeweils betrachteten Experiments ($j$).

Ergebnisse


Für dieses Anwendungsbeispiel wurden $106$ Startwert-Kombinationen innerhalb der oben angegebenen Modell-Parameter-Grenzen nach dem Saltelli-Sampling-Schema erzeugt und mittels Zielfunktion $J\left(\theta\right)$ (Gl. \ref{eq:costfun}) bewertet. Die nachfolgende Parameter-Identifikation mittels Minimierung der Zielfunktion ausgehend von geeigneten Startwert-Kombinationen konvergierte abschließend zu den in nachfolgender Tabelle aufgeführten Parameter-Beträgen.

Symbol Bezeichnung Betrag Einheit
$\mu_{\max}$ maximale spezifische Wachstumsrate $0.32$ $\text{h}^{-1}$
$\rho_{\max}$ maximale spezifische Produktbildungssrate $0.020$ $\text{h}^{-1}$
$K_{\text{S}}$ Halbsättigungskonstante Substrat $2.3$ $\text{g}_{\text{S}} \cdot \text{L}^{-1}$
$K_{\text{P}}$ Halbsättigungskonstante Produkt $1.6$ $\text{g}_{\text{S}} \cdot \text{L}^{-1}$
$K_{\text{I,S}}$ Inhibitionskonstante Substrat $9.1$ $\text{g}_{\text{S}} \cdot \text{L}^{-1}$
$K_{\text{I,P}}$ Inhibitionskonstante Produkt $0.10$ $\text{g}_{\text{P}} \cdot \text{L}^{-1}$
$k_{\text{d}}$ spezifische Biomasse-Zerfallsrate $0.024$ $\text{h}^{-1}$
$m$ spezifische Substrat-Maintenance-Rate $0.034$ $\text{g}_{\text{S}} \cdot \text{g}_{\text{X}}^{-1} \cdot \text{h}^{-1}$
$Y_{\text{X,S}}$ Ertragskoeffizient Biomasse/Substrat $0.60$ $\text{g}_{\text{X}} \cdot \text{g}_{\text{S}}^{-1}$
$Y_{\text{P,S}}$ Ertragskoeffizient Produkt/Substrat $0.053$ $\text{g}_{\text{P}} \cdot \text{g}_{\text{S}}^{-1}$

Nachfolgende Abbildung visualisiert, dass die experimentelle Datenbasis durch das parametrierte Bioprozess-Modell in allen wesentlichen Systemeigenschaften abgebildet werden kann (zum Erweitern anklicken):

Vergleich zwischen Messung und Simulation

Das parametrierte Bioprozess-Modell kann im nachfolgenden zur Analyse (z. B. Sensitivitätsanalyse) oder zur Optimierung (z. B. Ausbeute-Maximierung) des bestehenden Bioprozesses eingesetzt werden.

Referenzen


  1. [1] L. Gao, Y. Ren, Y. Ma, J. Lin, J. Lin: Modeling and simulation of production of metallothionein and red fluorescent fusion protein by recombinant Escherichia coli using graphical programming. In: Modeling, Programming and Simulations Using LabVIEW Software., IntechOpen, London, United Kingdom, 2011. doi: 10.5772/14091
  2. [2] D. Voet, J. G. Voet: Biochemistry. John Wiley & Sons, New York, 1990. ISBN: 0-471-61769-5.
  3. [3] A. Kremling: Grundlagen der mathematischen Modellierung. In: Kompendium Systembiologie, Vieweg+Teubner Verlag, 2012. doi: 10.1007/978-3-8348-8607-1_3.
  4. [4] J. A. Nelder, R. Mead: A simplex method for function minimization. In: The Computer Journal, 7(4):308–313, 1965. doi: 10.1093/comjnl/7.4.308
  5. [5] A. Raue, C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmüller, J. Timmer: Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. In: Bioinformatics, 25(15):1923–1929, 2009. doi: 10.1093/bioinformatics/btp358
  6. [6] B. Tang: Orthogonal array-based latin hypercubes. In: Journal of the American Statistical Association, 88(424):1392–1397, 1993. doi: 10.2307/2291282
  7. [7] I. M. Sobol: Distribution of points in a cube and approximate evaluation of integrals. Zh. Vych. Mat. Mat. Fiz., 7:784–802, 1967.
  8. [8] A. Saltelli: Making best use of model evaluations to compute sensitivity indices. In: Computer Physics Communications, 145:280–297, 1967. doi: 10.1016/S0010-4655(02)00280-1