Controladores biomoleculares

¿Qué mecanismo permite que la expresión génica sea tan robusta frente a perturbaciones externas y al ruido?

La expresión génica en células es notablemente robusta: es capaz de mantenerse dentro de rangos funcionales en todo momento a pesar de las fluctuaciones en las condiciones intracelulares, extracelulares y el ruido. Las células parecen regular la expresión de sus genes de manera que ''amortiguan'' estas fluctuaciones, lo que les permite mantener un comportamiento deseado con precisión y en el momento adecuado. Muchas vías de señalización exhiben este nivel de control, incluyendo la quimiotaxis (Yi et al., 2000), la respuesta osmótica (Muzzey et al., 2009) y varios procesos electrofisiológicos (De Palo et al., 2013), los cuales tienden a estabilizar su salida en torno a un valor objetivo específico incluso en presencia de perturbaciones persistentes. Entonces, ¿qué tipo de estrategia utilizan las células para lograr una expresión génica tan robusta?

Este tipo de regulación se logra gracias al control por feedback. En un esquema general de control (ver Figura 1, izquierda), este tipo de controlares generan una señal de control ($U$) basada en una métrica de error estimada al comparar la salida actual del proceso que se desea controlar ($Y$) con un valor objetivo deseado ($R$), también llamado referencia. El objetivo principal del controlador es minimizar este error. De este modo, si ocurre una perturbación en el proceso o si el entorno es ruidoso, los niveles de expresión génica son corregidos de vuelta hacia su valor deseado (ver Figura 1, derecha).

En esta clase exploraremos cómo diseñar uno de los controladores de feedback más simples basado en una reacción de secuestramiento molecular: el controlador proporcional. Utilizaremos el principio de separación de escalas de tiempo y análisis en estado estacionario para entender cómo podemos reducir el error entre el nivel deseado de expresión génica y el nivel real. También realizaremos un análisis de estabilidad local para derivar condiciones matemáticas que garanticen asintóticamente que el controlador permanezca acotado. Finalmente, examinaremos los principios de diseño para construir un controlador de feedback basado en una reacción química que exhiba ultrasensibilidad.

Fundamentos del diseño de un controlador proporcional

El controlador proporcional es el tipo más simple de controlador de feedback, en el cual la señal de control es directamente proporcional al error entre la salida y la reference. Como ejemplo, diseñaremos un controlador basado en secuestramiento molecular, donde una especie $R$ actúa como referencia y produce $U_1$, que interactúa con $U_2$ mediante una reacción de secuestramiento. En cuanto al proceso (a veces llamado también planta), consideraremos un sistema de expresión en cascada de dos especies, donde una especie química $W$ produce una especie $Y$, la cual es luego detectada por el controlador para formar el sistema en lazo cerrado mostrado en la Figura 2 (izquierda).

Comencemos escribiendo las reacciones químicas que describen al sistema de lazo cerrado. Cabe resaltar que estamos ignorando la dinámica del complejo $C$ mostrado en la Figura 2 debido a que no afecta la dinámica del sistema completo en el régimen de secuestramiento rápido (ver Lecture Notes #1).

\begin{align} R &\xrightarrow{\alpha_1} U_1 && U_1 \xrightarrow{\delta} \emptyset && \text{Producción y degradación of $U_1$}\\ Y &\xrightarrow{\alpha_1} U_2 && U_2 \xrightarrow{\delta} \emptyset && \text{Producción y degradación of $U_2$}\\ U_1 + U_2 &\xrightarrow{\gamma} \emptyset && && \text{Secuestramiento molecular} \\ U_1 &\xrightarrow{\theta_1} U_1 + W && W \xrightarrow{\delta} \emptyset && \text{Señal de control y degradación} \\ W &\xrightarrow{\theta_2} W + Y && Y \xrightarrow{\delta} \emptyset && \text{Producción y degradación of $Y$} \end{align}

A continuación, modelaremos estas reacciones químicas utilizando la ley de acción de masas. Sin embargo, vemos que existe un número considerable de parámetros cinéticos ($\alpha_1$, $\theta_1$, $\theta_2$, $\delta$, $\gamma$), lo cual puede dificultar nuestro análisis matemático. Por ello, utilizaremos una versión adimensional de las ecuaciones diferenciales para reducir la cantidad de parámetros involucrados en el análisis (ver Technical Note #4). Como resultado, trabajaremos con un sistema un poco más sencillo:

\begin{eqnarray} \frac{d}{d\tau} u_1 &=& \alpha_1 r - u_1 - \frac{1}{\xi} u_1 u_2 \\ \frac{d}{d\tau} u_2 &=& \alpha_1 y - u_2 - \frac{1}{\xi} u_1 u_2 \\ \frac{d}{d\tau} w &=& \theta_1 u_1 - w \\ \frac{d}{d\tau} y &=& \theta_2 w - y \end{eqnarray}

Definamos el error en estado estacionario ($e$) como la diferencia entre la referencia, dada por la concentración (constante) de la especie $R$, y la especie $Y$, de modo que $e = r - \bar y$. Ahora exploremos el régimen cinético en el que la tasa de secuestro es extremadamente alta, es decir, cuando $\xi \to 0$. Al resolver el estadio estacionario de las ecuaciones en función del error, obtendremos la siguiente expresión para $e$:

\begin{eqnarray} \lim_{\xi \to 0} e = \frac{r}{1 + \alpha_1 \theta_1 \theta_2} \end{eqnarray}

En otras palabras, si aumentamos la ganancia del controlador ($\alpha_1$), podremos reducir el error en estado estacionario. Sin embargo, es importante preguntarnos: ¿cuánto podemos aumentar esta ganancia sin desestabilizar el proceso? Como un primer experimento, realicemos algunas simulaciones para observar cómo responde el sistema a medida que incrementamos $\alpha_1$, buscando el valor más alto que aún mantenga estabilidad.

Mediante prueba y error, encontramos que cuando $\alpha_1$ se aproxima o supera el valor de 8, el sistema comienza a oscilar. Esto indica que el proceso ha pasado a un régimen dinámico distinto y ya no converge al valor de la referencia. Esto es un ejemplo de la compensación entre reducir el error en estado estacionario y mantener la estabilidad del sistema. Para formalizar esta observación, podemos utilizar un análisis de estabilidad local para derivar una condición que garantice estabilidad asintótica. Al operar en el régimen donde el secuestramiento es muy rápido ($\xi \to 0$), el análisis de los autovalores del polinomio característico (ver Technical Note #3 y Lecture Notes #2) nos lleva a la condición:

\begin{eqnarray} 0 < \alpha_1 \theta_1 \theta_2 < 8 \end{eqnarray}

, que es precisamente lo que encontramos con nuestras simulaciones!

Reduciendo el error a cero: diseñando un controlador (casi) integral

En 2016, Briat et al. formalizaron el mecanismo mediante el cual ciertos sistemas biológicos pueden lograr adaptación perfecta; es decir, la capacidad de ''recuperarse'' de una perturbación y volver exactamente a su estado original, sin dejar de responder a cambios posteriores. Utilizando conceptos de teoría de control, los autores propusieron el controlador integral antitético, un circuito basado en la reacción de secuestramiento molecular. Demostraron tanto de manera matemática (Briat et al., 2016a) como experimental (Briat et al., 2016b) que esta circuito puede generar un error en estado estacionario igual a cero, siempre que la degradación de las especies sea despreciable (~ 0).

Dado que esta acción integral depende de la tasa de degradación, debemos adaptar nuestras ecuaciones anteriores para reflejar que las especies del controlador ($U_1$ y $U_2$) se degradan de manera distinta a las especies del proceso ($W$ y $Y$). Para representar esto, introducimos un parámetro adimensional $\zeta$ que multiplica los términos de degradación de $u_1$ y $u_2$, y exploramos el régimen cinético en el que $\zeta \to 0$.

\begin{eqnarray} \frac{d}{d\tau} u_1 &=& \alpha_1 r - \zeta u_1 - \frac{1}{\xi} u_1 u_2 \\ \frac{d}{d\tau} u_2 &=& \alpha_1 y - \zeta u_2 - \frac{1}{\xi} u_1 u_2 \\ \frac{d}{d\tau} w &=& \theta_1 u_1 - w \\ \frac{d}{d\tau} y &=& \theta_2 w - y \end{eqnarray}

Nota: esta operación es equivalente a asumir degradaciones distintas (a saber, $\delta$ para el controlador y $\phi$ para el proceso), y repetir los pasos de adimensionalización con respecto a $\phi$.

Ahora que hemos confirmado las condiciones bajo las cuales el controlador integral antitético funciona correctamente, estudiaremos más de cerca dos propiedades importantes: el tracking de referencia, que es la capacidad del sistema para ajustarse a un nuevo valor de referencia, y la robustes ante perturbaciones, que es su capacidad para volver al estado anterior después de haber sido perturbado.

Control adaptativo

Al diseñar controladores de feedback, uno de los desafíos más importantes es que generalmente necesitamos entender cómo se comporta el sistema. En la práctica, esto significa utilizar ecuaciones que describen el proceso que queremos controlar (a esto también se le denomina tener un entendimiento mecanístico del sistema). Sin embargo, la mayoría de las vías de señalización en biología aún no se entienden por completo. Por ello, puede resultar difícil elegir una buena señal de referencia o ajustar correctamente los parámetros del controlador, lo que lleva a ciclos iterativos de prueba y error. Esto nos plantea una pregunta: ¿es posible construir un controlador de feedback incluso cuando solo conocemos parcialmente el comportamiento del sistema? Los controladores adaptativos son una forma de hacerlo. En lugar de depender de una referencia externa, este tipo de controladores son capaces de estimarla a partir de la propia dinámica del sistema, ajustándose a medida que el sistema cambia (ver Figura 3).

Diseño de un controlador adaptativo biomolecular usando el motivo incoherent feedforward loop (IFFL)

Una de las estrategias más directas para implementar control adaptativo se basa en la acción derivativa. En 2004, Pyragas et al. implementaron un controlador adaptativo (electrónico) basado en un filtro pasa-bajas y un operador derivativo. Esta estrategia no afecta el equilibrio del sistema, ya que, si la salida del controlador es proporcional a la derivada del proceso, cuanto más cerca se encuentra del estado estacionario, menor será la acción del controlador. Una vez que el sistema alcanza su estado estacionario, la acción del controlador desaparece, ya que la derivada converge a cero (Pyragas et al., 2004).

Ahora estudiaremos cómo implementar un diseño similar utilizando reacciones químicas. Para lograr un funcionamiento análogo, consideramos el motivo conocido como *incoherent feedforward loop* (IFFL), en el que una entrada $Y$ produce simultáneamente las especies $U_1$ y $U_2$, y donde $U_2$ inhibe la acción de $U_1$. Consideremos ahora una reacción de secuestro molecular entre las especies $U_1$ y $U_2$. Supongamos además que otra especie intermedia $X$ produce $U_1$, y que tanto $X$ como $U_2$ son producidas por la entrada $Y$ (ver Figura x). Esta estrategia implementa esencialmente una topología del tipo IFFL, la cual podemos modelar mediante el siguiente sistema de ecuaciones diferenciales ordinarias:

\begin{eqnarray} \frac{du_1}{dt} &=& kx - \delta u_1 - \gamma u_1 u_2 \\ \frac{du_2}{dt} &=& \rho y - \delta u_2 - \gamma u_1 u_2 \\ \frac{dx}{dt} &=& kx - \delta x \end{eqnarray}

Para nuestro análisis, consideremos el complejo $C$ que se forma cuando $U_1$ y $U_2$ interactúan mediante secuestramiento, pero que usualmente ignoramos ya que no juega un rol en la dinámica del sistema de lazo cerrado cuando el secuestramiento es muy rápido. Podemos modelar la dinámica de $C$ con la siguiente ecuación diferencial:

\begin{eqnarray} \frac{dc}{dt} &=& \gamma u_1 u_2 - \delta c \end{eqnarray}

Ahora, consideremos las ecuaciones de conservación de masa para $u_1$ y $u_2$. Sea $u_1^T$ la concentración total de $u_1$ (ya sea libre o unida a $U_2$ formando el complejo), de modo que $u_1^T = u_1 + c$. Lo mismo para la especie $U_2$: sea $u_2^T = u_2 + c$. Si reescribimos nuestras ecuaciones diferenciales anteriores en términos de $u_1^T$ y $u_2^T$, notaremos que son lineales!

\begin{eqnarray} \frac{d u_1^T}{dt} &=& kx - \delta u_1^T \\ \frac{d u_2^T}{dt} &=& \rho y - \delta u_2^T \\ \frac{d x}{dt} &=& \theta y - \delta x \end{eqnarray}

De nuestro Lecture Note #1, recordamos que podemos escribir la dinámica de $u_1(t)$ como $\max (0, u_1^T - u_2^T)$ en el régimen de secuestramiento rápido. Podemos reescribir esta expresión como $u_1(t) \approx \max \{ 0, \mathcal{L}^{-1} [U_1^T(s) - U_2^T(s)] \}$, donde $\mathcal{L}^{-1}$ es la transformada inversa de Laplace, y $U_1^T(s)$ y $U_2^T(s)$ son las transformadas de Laplace de $u_1^T(t)$ y $u_2^T(t)$, respectivamente. Ahora, al aplicar la transformada de Laplace a nuestras ecuaciones diferenciales anteriores y reorganizar los términos, encontramos que

\begin{eqnarray} u_1(t) \approx \max \left \{ 0, \mathcal{L}^{-1} \left ( \frac{k \theta}{(s+\delta)^2} \left [1 - \frac{\rho \delta}{k \theta} - \frac{\rho}{k \theta}s \right ] Y(s) \right) \right \} \end{eqnarray}

A partir de esta expresión para la dinámica aproximada de $u_1(t)$, definamos la métrica $r = \rho \delta / k \theta$. Cuando $r = 1$, notamos que el término $1 - r$ desaparece de la ecuación. Por lo tanto, nos queda $k \theta / (s + \delta)^2$, que representa un filtro pasa-bajas, y $- \rho/k \theta \, s$, que representa un operador derivativo que responde a gradientes negativos; exactamente los dos operadores que necesitábamos para el controlador adaptativo! Por esta razón, llamaremos a $r$ la métrica adaptativa ideal y exploraremos sus propiedades mediante algunas simulaciones.

Desestabilización del estado estacionario de un gen auto-inhibidor

Como ejemplo, consideremos un gen que expresa un represor transcripcional que inhibe su propia transcripción. Si asumimos una dinámica sin ningún retraso adicional, tendremos un modelo simple que representa un sistema biológico que ha perdido su comportamiento oscilatorio debido a una enfermedad, como ocurre con la disrupción de los ritmos circadianos en el cáncer (Zambrano et al., 2016) y en enfermedades metabólicas (Gabriel et al., 2021). Volvamos a adimensionalizar el sistema, lo que nos lleva al siguiente conjunto de ecuaciones diferenciales para el gen que inhibe su propia transcripción conectado al controlador adaptativo:

\begin{eqnarray} \frac{dy}{dt} &=& \alpha \frac{1}{1 + y^m} - y + \beta u_1 \\ \frac{du_1}{dt} &=& kx - u_1 - \frac{1}{\xi} u_1 u_2 \\ \frac{du_2}{dt} &=& \rho \frac{1}{1 + y^m} - u_2 - \frac{1}{\xi} u_1 u_2 \\ \frac{dx}{dt} &=& kx - x \end{eqnarray}

donde el término $\beta u_1$ representa la señal de control del controlador adaptativo.

Dado que el controlador adaptativo es capaz de generar oscilaciones sostenidas, surge una pregunta natural: ¿son esas oscilaciones robustas frente a perturbaciones? O, de forma aún más interesante, ¿son las propiedades dinámicas de esas oscilaciones, como su frecuencia y amplitud, robustas frente a perturbaciones? Antes de explorar estas preguntas mediante simulaciones, definamos el período como el inverso de la frecuencia, calculado como el intervalo promedio entre picos consecutivos. De manera similar, definamos la amplitud como la diferencia promedio entre los valores de esos picos.

Lecturas recomendadas

[1] Del Vecchio, D., Dy, A.J. and Qian, Y. (2016) ‘Control theory meets synthetic biology’, Journal of The Royal Society Interface, 13(120), p. 20160380. doi:10.1098/rsif.2016.0380.

[2] Cuba Samaniego, C. and Franco, E. (2021) ‘Ultrasensitive molecular controllers for quasi-integral feedback’, Cell Systems, 12(3). doi:10.1016/j.cels.2021.01.001.

[3] Filo, M., Kumar, S. and Khammash, M. (2022) ‘A hierarchy of biomolecular proportional-integral-derivative feedback controllers for robust perfect adaptation and dynamic performance’, Nature Communications, 13(1). doi:10.1038/s41467-022-29640-7.

[4] Britto Bisso, F. et al. (2024) Design of a biomolecular adaptive controller to restore sustained periodic behavior [Preprint]. doi:10.1101/2024.10.11.617962.