Dinámica de un Quimiostato no autónomo con retardo
Full text
DINÁMICA DE UN QUIMIOSTATO NO AUTÓNOMO CON RETARDO Trabajo fin de máster presentado por: Javier López de la Cruz bajo la dirección de: D. Tomás Caraballo Garrido Dpto. de Ecuaciones Diferenciales y Análisis Numérico VoBodel director Fdo.: Tomás Caraballo Garrido. MÁSTER UNIVERSITARIO EN MATEMÁTICA AVANZADA FACULTAD DE MATEMÁTICAS UNIVERSIDAD DE SEVILLA
A mis padres y a mi hermana, porque gracias a ellos soy quien soy, por todo el ánimo y apoyo recibido, por toda la confianza que ponen en mi, por ser el pilar base de mi vida. A Tomás, mi director, porque gracias a él esto ha sido posible, por todo lo que, día a día, me enseña, por todas las oportunidades que me brinda, por abrirme hueco en este apasionante mundo. A Julia Clara, mi gran amiga, por todos los ratos compartidos, por ser mi compañera en este viaje, por entrar a formar parte de mi vida, porque sin ella nada habría sido igual.
Agradecimientos Quiero aprovechar estas primeras páginas que dan comienzo a este trabajo para hacer partícipes de él a algunas de las personas que, junto a mí, han sido testigos y me han ayudado, tanto en el terreno académico como en el personal, a que todo este proyecto, que he hecho con una gran ilusión, pueda llegar a ver la luz por fin. Merecen mención especial mis padres, Toñi y Longi, por haberme dado la vida y con ello la oportunidad de disfrutar haciendo cada día lo que hago, lo que me apasiona, aquello que deseo seguir realizando el resto de mi vida. Ellos siempre están ahí, a pesar de los trescientos kilómetros de distancia que nos separan cada vez que vuelvo a Sevilla “a formarme para intentar ayudar al mundo con, al menos, pequeños granitos de arena”, tal y como les digo tantas veces cuando me despiden en la estación. Siempre me han apoyado, me han animado a hacer lo que me gusta, aunque eso supusiera no tenerme en casa como hasta el momento y poder estar conmigo cada día, y siempre han sabido proporcionarme la fuerza necesaria, tanto en los momentos felices como en aquéllos que no lo eran tanto. Y lo más importante de todo: si hoy día soy yo, es sin duda gracias a ellos. También merece mención especial mi hermana pequeña, Ana Belén, o mi Anita, como yo le digo, que aguanta mis frases hablando con términos matemáticos, se ríe y las repite, sacándome algunas de las sonrisas más sinceras. Porque su fortaleza me sorprende cada día, porque siempre me apoya en todo lo que decido hacer y porque siempre habla con mucho orgullo de su hermano matemático que estudia en Sevilla. Sin duda, mis padres y ella, quienes más tiempo tienen que soportar todas mis explicaciones sobre lo que voy aprendiendo en la facultad, son mis referentes y juntos forman el pilar base de mi vida. Otra persona muy especial para mi es Tomás, mi director, quien me enseña cada día y me transmite todo lo que sabe de una forma que admiro; todos los que me conocen saben que es para mi un modelo a seguir, no sólo por la brillantez y el entusiasmo en sus explicaciones, sino también en lo que respecta a la vida cotidiana. No tendré nunca tiempo ni palabras suficientes para agradecerle como se merece todo lo que ha hecho y hace por mi, pues todas las oportu- i
ii AGRADECIMIENTOS nidades que se me brindan en este campo que me fascina son gracias a él. Espero poder llegar algún día a hacer lo que él hace y a tener alumnos que me admiren tanto como yo a él. Pero no sólo eso, y es que no puedo acabar sin agradecerle su disponibilidad cada vez que necesito su ayuda ¡y eso que algunas veces puedo llegar a resultar algo pesado! Ahora me gustaría hacer un breve recorrido por todas las personas que han estado junto a mi en los últimos seis años de mi estancia en Sevilla, a todos mis amigos y trabajadores de la Residencia Universitaria Armendáriz, donde nací como universitario; allí me acogieron como a un hijo tras haber salido del cascarón por primer vez. No puedo pasar sin mencionar a uno de los trabajadores, a un gran amigo que siempre ha estado ahí para ayudar a todo el que lo necesitaba porque estaba pasando por un mal período de su vida, a nuestro Ramón. Allí pasé los tres primeros años de mi vida universitaria; sin embargo, los tres últimos los he pasado en la Residencia Universitaria Rector Estanislao del Campo, donde he conocido a muchísimos amigos también, personas estupendas y trabajadores extraordinarios. Tengo que agradecer a todas las personas con las que he vivido estos seis años el hecho de haber creado un clima universitario, lleno de experiencias inolvidables, que siempre permanecerá en nuestros corazones. Especialmente quiero mencionar a Julia Clara, con quien he compartido muy de cerca estos dos últimos años, con quien he pasado tantos ratos de estudio; ella ha sido una de las personas que ha vivido más de cerca mis horas de dedicación y ha estado dispuesta a ayudarme en mis ratos menos alegres. Además, es una de las personas que más se ha interesado por saber qué era lo que estaba haciendo, junto con mis padres y mi hermana, así como conocer el transcurso de los progresos. Estas líneas son para ella, para agradecerle todos los ratos que hemos pasado juntos y todos los abrazos que me ha dado, en definitiva, un agradecimiento a mi mejor amiga. Finalmente, quiero dar mi agradecimiento a mis abuelos, por el ánimo que siempre me dan y a todos los profesores que he tenido a lo largo de estos años, especialmente a los profesores del Departamenteo de Ecuaciones Diferenciales y Análisis Numérico, por su amabilidad y simpatía cada vez que me ven. Javier López de la Cruz.
Prólogo Un quimiostato es un mecanismo de laboratorio que consiste en tres recipientes interconectados y es usado para el crecimiento de microorganismos en un ambiente de cultivo. En su forma más básica, el excedente del primer recipiente pasa al segundo y el excedente del segundo pasa al tercero. El primer recipiente se denomina botella de alimentación y contiene todos los nutrientes requeridos para el crecimiento de los microorganismos. Suponemos que todos los nutrientes se aportan de forma abundante excepto uno de ellos, que denominaremos nutriente limitante. El contenido de la botella de alimentación es bombeado al segundo recipiente, al que denominamos recipiente de cultivo, con un ratio constante. Los microorganismos se alimentan de los nutrientes provenientes de la botella de alimentación y crecen en el recipiente de cultivo. El recipiente de cultivo se remueve de forma continua para que todos los microorganismos tengan la misma probabilidad de acceder a los nutrientes. El contenido del recipiente de cultivo es bombeado entonces al tercer recipiente, al que denominamos recipiente colector. Naturalmente tal recipiente contiene nutrientes, microorganismos y los productos generados por tales microorganismos [1]. Figura 1: Representación gráfica de un quimiostato. iii
iv PRÓLOGO El quimiostato juega un papel importante en estudios ecológicos [7, 8, 9, 10, 11, 12, 13, 14]. Con algunas modificaciones, también es usado como modelo para procesos de tratamiento de aguas residuales [15, 16]. El modelo del quimiostato puede considerarse como el punto de partida para muchas variantes que cedan el paso a modelos biológicos más realistas, por ejemplo, problemas de microorganismos genéticamente alterados [17, 18] y modelos de antibióticos [37]. Más literaturas en cuanto a variantes y análisis de modelos del tipo del quimiostato pueden encontrarse en [19, 20, 21] y sus referencias. En el modelo más simple de quimiostato, la disponibilidad de nutriente y su ratio de administración se suponen fijados. Sin embargo, la disponibilidad de un nutriente en un sistema natural normalmente depende del ratio de consumo del nutriente y de la concentración del mismo, lo que conduce a un sistema dinámico no autónomo. Otro supuesto básico en el modelo más simple de quimiostato es que el ratio de flujo sea suficientemente rápido como para no permitir el crecimiento en las paredes del recipiente. Sin embargo el crecimiento en la pared tiene lugar cuando el ratio de limpieza no es suficientemente rápido, lo que da lugar a los problemas de bioreactores. Los estudios de modelos de quimiostato tratados como sistemas dinámicos no autónomos están mucho más limitados hasta la fecha, por ejemplo, Smith y Thieme introdujeron persistencia práctica para sistemas dinámicos no autónomos con el modelo más simple de quimiostato, ejemplo en [22] cuando el ratio de limpieza depende del tiempo. Antes de dar un breve resumen sobre lo que trataremos en este trabajo, hablaremos sobre la analogía entre un quimiostato y un lago, de hecho una forma sencilla de entenderlos para los biólogos es mirarlos como si de un simple lago se tratasen. De esta forma, los modelos de quimiostatos son usados extensamente para representar el crecimiento de especies en un lago donde los organismos, como por ejemplo algas, se alimentan de nutrientes con crecimiento limitado, como pueden ser el nitrógeno y el fósforo. La analogía entre un lago y un quimiostato puede verse de forma clara en la tabla que aparece más adelante. La disponibilidad de un nutriente en un sistema natural como es un lago depende de la afluencia de nutriente. Las comunidades de algas, por ejemplo, sobreviven incluso con niveles muy bajos (casi indetectables) de nutrientes, en contra de la opinión de que perecen debido a insuficiencias. Pero hay un crecimiento oscilatorio y bajo. Para representar este fenómeno de crecimiento oscilatorio en las ecuaciones, los investigadores han intentado varios caminos.
PRÓLOGO v Propiedad Lago Quimiostato Nutriente Las especies en un lago reciben los nutrientes a través de las corrientes que fluyen en o por la regeneración durante la primavera o el otoño. El nutriente se suministra a través de una entrada. Muerte Las especies se extinguen, ya que continuamente caen fuera de las capas luminosas a las capas inferiores de la columna de agua. El nutriente y la especie se lavan fuera del sistema a través de una salida. Escasez de nutriente Durante el verano no hay suficiente aportación de nutriente y debido a la sequía del lago, algunos nutrientes tales como fósforo, nitrógeno o vitamina B12 están menos disponibles. La aportación de uno o más nutrientes (esencial) está controlada por el experimentador. Cuadro 1: Tabla de comparación entre lago y quimiostato. En el primero de todos, se introducen una cantidad de nutriente y un ratio de limpieza [23, 24, 25, 26] variables en el modelo, los cuales se asumen fijos en un principio, debido a que la disponibilidad de los nutrientes (I) y también su ratio de aportación (D) en un lago dependen de la estación del año en la que nos encontremos. Para hacernos una idea, los investigadores han hecho variar DeIde forma periódica de acuerdo con el tiempo. Esto puede explicar el crecimiento oscilatorio (coexistencia en caso de competiciones). Pero asumir periodicidad es muy sencillo, ya que la influencia puede variar continuamente durante una estación en particular o puede permanecer regulada. En otras palabras, los parámetros IyDpueden variar sobre un ciclo o pueden estar fijados (constantes) o con variaciones marginales durante una parte específica del ciclo. Por tanto, su influencia en el sistema debe entenderse para estaciones específicas o algunos períodos de tiempo. Otra posibilidad que los matemáticos manejan consiste en intentar introducir otros niveles tróficos. Así, el modelo se modifica para incluir un nutriente, un microorganismo que se alimenta de tal nutriente, y dos competidores que se alimentan del microorganimo [27, 28, 29, 30, 31, 32, 33, 34]. También se demuestra que hay coexistencia de todas las especies siempre y cuando se establezca la presencia de unos ciclos límite bajo condiciones apropiadas sobre los parámetros, pero no es posible introducir un depredador (u otro nivel trófico) en un sistema cerrado tal como un quimiostato cuando estabamos estudiadno el crecimiento de un nutriente - un microorganismo (o competidor) sólo.
Capítulo 1 Modelo con un microorganismo. 1.1. Introducción. Iniciamos este trabajo recalcando la relación clásica que tiene lugar entre las matemáticas y la biología, antes de comenzar con la explicación, el estudio y el análisis detallado del modelo que nos ocupa. Ponemos de manifiesto por tanto que la aplicación de las matemáticas a los problemas cuyo origen reside en la biología no son nuevos, recordemos los modelos Malthusianos y de Lotka-Volterra, entre otros. El modelo, x′(t) = rx(t), donde ′=d/dt yres un coeficiente que refleja el ratio de crecimiento, describe el crecimiento exponencial de una población xde individuos en cualquier instante de tiempo t. La introducción de una constante de capacidad de carga del medio, como bien pueden ser la comida, la luz, el espacio, etcétera, que toman control sobre el crecimiento de las especies, modifica el modelo anteriormente citado logrando llegar a un modelo más realista de la forma x′=rx 1−x k, donde kes dicha constante de capacidad de carga del medio. Un modelo del tipo Lotka-Volterra que describe la competición entre dos especies xey(depredadores) por la misma comida (presa) viene dado por x′=rx 1−x k1, y′=ry 1−y k2.(1.1) 5
6CAPÍTULO 1. MODELO CON UN MICROORGANISMO. Nótese que este modelo es un modelo sin interferencia. Podemos tambier considerar un modelo con interferencia, como el que sigue x′=rx 1−x k1−λ1y, y′=ry 1−y k2−λ2x.(1.2) En los modelos (1.1) y (1.2), k1,k2son constantes de capacidad de carga y λ1,λ2son constantes de interferencia, todas ellas constantes positivas. Es posible medir los coeficientes de interferencia λ1yλ2que deciden el resultado de la competición sólo cuando los organismos crecen juntos, es decir, mientras el experimento está siendo realizado. De esta forma se descarta la posibilidad de predecir el resultado del experimento antes de su realización. Este es un inconveniente del modelo de Lotka-Volterra. Además, no hay representación para el crecimiento de la comida (fuente) gracias a la cual la población crece, debido a que se alimenta de ella. La inclusión de la ecuación que describe el crecimiento de la fuente (presa) nos permite predecir los cambios correspondientes en el crecimiento de la población consumidora (depredaror). La importancia de los modelos de quimiostatos reside en la inclusión de la ecuación que describe el crecimiento de la fuente junto con las ecuaciones que describen la competición entre los organismos que intervienen en la misma. Esta observación ha incitado a los investigadores a estudiar las ecuaciones de los modelos de quimiostatos de forma intensiva. De esta forma, analizando las ecuaciones del modelo, el resultado de la competición que tiene lugar en un quimiostato se puede predecir y dichas predicciones pueden ser confirmadas por el experimento. Llegados a este punto, comenzaremos nuestro camino hacia la construcción del modelo matemático. Para ello, explicaremos a continuación los principios biológicos bajo los que el modelo será formulado. 1.2. Principios biológicos. El modelo con el que trabajaremos a lo largo de toda esta memoria consiste en un microorganismo que se alimenta de un único nutriente cuyo crecimiento es limitado. Denotaremos por xal nutriente cuyo crecimiento es limitado y por y, al microorganismo que se alimenta del nutriente x. A continuación enunciaremos algunos supuestos que nos serán de ayuda a la hora de formular nuestro modelo. ⋆Todos los nutrientes, excepto x, se encuentran en total disponibilidad, pero no suplen la posible ausencia de x.
1.2. PRINCIPIOS BIOLÓGICOS. 7 ⋆Contamos con una fuente externa desde la que el microorganismo recibe el nutriente x junto con otros nutrientes. ⋆El nutriente es administrado con un ratio y el nutriente aportado tiene una concentración en el medio de crecimiento. ⋆El nutriente es uniformemente distribuído en el medio de crecimiento, de tal forma que todos los microorganimos tienen igual acceso al nutriente disponible, esto es, no hay disparidad en la distribución. ⋆Hay un excedente en el sistema donde el material que se encuentra en exceso, ya sean nutrientes, o bien microorganismos, y algunos otros productos, es eliminado continuamente. Si paramos a pensar un poco, el último supuesto citado es razonable ya que en cualquier sistema tal como un dispositivo de laboratorio, el volumen es fijado. De esta forma, debemos contar con un excedente al permitir disponibilidad de material de forma continua. Además, para preservar el volumen del medio de crecimiento (fijado), debemos asumir que el ratio de aporte y el ratio de eliminación de los contenidos del medio de crecimiento son el mismo. A continuación enumeramos una serie de aspectos a tener en cuenta. 1. Los microorganismos consumen nutriente continuamente y el ratio de dicho consumo es una constante. 2. El consumo posee una saturación, en otras palabras, la suplencia ilimitada de nutrientes no implica el consumo ilimitado de nutriente por las especies. 3. El crecimiento de los microorganismos es proporcional a su consumo. 4. Todos los factores externos tales como temperatura, presión, etc. son apropiados para el crecimiento y no afectan al sistema. De esta forma, tenemos algunos parámetros reales y relaciones funcionales en el sistema. Veamos una enumeración de los mismos. ⋆I, cantidad de nutriente disponible en el sistema en cualquier instante de tiempo que es una constante positiva. ⋆D, ratio según el que el nutriente es aportado y también ratio según el que los contenidos del medio de crecimiento son eliminados. También es una constante positiva. ⋆a, ratio de consumo máximo de nutrientes y, por tanto, ratio de crecimiento de los microorganismos, constante positiva.
8CAPÍTULO 1. MODELO CON UN MICROORGANISMO. ⋆U, respuesta funcional de los microorganismos que describe cómo se produce el consumo de nutrientes. Ya estamos en condiciones de introducir nuestro modelo matemático que describe la dinámica de un sistema con consumo de nutriente limitado. 1.3. Formulación del modelo matemático. Denotamos por x(t),y(t)a las concentraciones de nutriente y microorganismos en un instante de tiempo tcualquiera, respectivamente. De esta forma, contamos en un principio con un término de ganancia: partiendo de I, se suministra al sistema, en cada instante de tiempo, un ratio Dde nutriente, por tanto, el ratio de incrementeo de la concentración de nutrientes en el instante tes DI. Además, tenemos los siguientes términos de pérdida: Una cantidad x(t)es eliminada del sistema con un ratio D. Además, el sistema pierde Dx(t)nutrientes cada instante de tiempo. Cada microorganismo yconsume U(x(t)) nutrientes. Además, U(x(t))y(t)es la concentración de nutrientes consumida por ycon un ratio a. De esta forma, la concentración de nutrientes viene dada por aU(x(t))y(t). Así, tenemos la siguiente ecuación: ratio de cambio de la concentración de nutrientes en el medio de crecimiento = =ratio de aportación de nutrientes −ratio de eliminación −ratio de consumo. Matemáticamente, dx(t) dt =DI −Dx(t)−aU(x(t))y(t).(1.3) Para el crecimiento de los microorganismos yasumimos que el consumo de nutrientes implica el crecimiento de y. Además, el ratio con el que tal nutriente es consumido, es decir, aU(x(t))y(t) es también el ratio con el que crece el microorganismo y. Por otro lado, el microorganismo es eliminado del sistema con un ratio constante D, con lo que Dy(t)denota la cantidad de microorganismos eliminados del sistema. Ahora tenemos la siguiente ecuación: ratio de cambio del microorganismo en cualquier instante de tiempo = =ratio de crecimiento del microorganismo −ratio de eliminación del microorganismo.
1.4. PROPIEDADES DE LAS SOLUCIONES. 9 Matemáticamente, dy(t) dt =aU(x(t))y(t)−Dy(t).(1.4) Las ecuaciones (1.3) y (1.4) describen, en su conjunto, el sistema dinámico de consumo con fuente limitada que intentamos conseguir en el quimiostato. Obtenemos entonces el siguiente sistema de ecuaciones diferenciales: dx(t) dt =DI −Dx(t)−aU(x(t))y(t), dy(t) dt =aU(x(t))y(t)−Dy(t).(1.5) 1.4. Propiedades de las soluciones. Debemos comenzar esta sección notando que las características del sistema (1.5) son esencialmente decididas por el término no lineal del miembro derecho, esto es, U(x).Use conoce en la literatura como la respuesta funcional o función de consumo. Hay algunos supuestos básicos sobre la función de consumo U: [0,∞)−→ [0,∞), que son los siguientes: 1. U(0) = 0,U(x)>0para x > 0. 2. l´ımx→∞ U(x) = L, con 0< L < ∞. 3. Ues continuamente diferenciable. 4. Ues monótona creciente. Los dos primeros supuestos son esenciales para cualquier función de consumo. El tercer supues- to proporciona una condición local de Lipschitz, lo que garantiza la existencial local de solución del sistema (1.5). El cuarto supuesto lo tendremos en cuenta, a pesar de no ser estrictamente necesario, ya que nos ayudará a concluir algunos de los resultados que se demuestran. Las siguientes funciones prototipo, que satisfacen los cuatro supuestos citados anteriormente, pueden tomarse como funciones de consumo a la hora de estudiar modelos biológicos tales como (1.5). Estas funciones son las que presentamos a continuación. U(x) = x: función de consumo Lotka-Volterra oHolling tipo-I. U(x) = x m+x,m > 0: función de consumo Michaelis-Menten oHolling tipo-II. U(x) = x2 (p+x)(q+x),p, q > 0: función de consumo senoidal oHolling tipo-III.
10 CAPÍTULO 1. MODELO CON UN MICROORGANISMO. La función de consumo U(x) = xno satisface el segundo supuesto, sin embargo, es usada con frecuencia en modelos del tipo Lotka-Volterra. Claramente los dos primeros supuestos garantizan la existencia de una constante positiva L > 0 tal que U(x)< L, para todo x∈[0,∞).(1.6) Además, el tercer supuesto junto con (1.6) garantizan la unicidad de solución del sistema (1.5); dichas soluciones son extensibles en su intervalo maximal de existencia para cualquier conjunto de condiciones iniciales. En lo que sigue, vamos a usar normalmente la función de consumo U(x) = x m+x, por lo que comprobaremos a continuación que, efectivamente, verifica los cuatro supuestos que ya hemos tratado. 1. Primer supuesto: U(0) = 0 m= 0.X 2. Segundo supuesto: l´ım x→∞ U(x) = l´ım x→∞ x m+x= 1 =: L∈(0,∞).X 3. Tercer supuesto: Debemos notar que puede parecer que la función de consumo posee un problema de discontinuidad en x=−m, sin embargo conviene recordar que mes una constante positiva, lo que indicaría una cantidad negativa de nutrientes, lo cual no es posible. Además, cualquier derivada de la función de consumo posee en el denominador una potencia de m+x, siendo x > 0ym > 0, con lo que se verifica este supuesto. X 4. Cuarto supuesto: U′(x) = m (m+x)2>0, de donde se deduce que Ues monótona creciente. X A continuación, vamos a demostrar que las soluciones del sistema (1.5) que comienzan con condiciones iniciales no negativas, se mantienen no negativas a lo largo del tiempo. Realmente este hecho es un requisitio que deben verificar xey, ya que están representando poblaciones,
1.4. PROPIEDADES DE LAS SOLUCIONES. 11 con lo que no pueden ser negativas en ningún caso. Teorema 1.4.1 Todas las soluciones del sistema (1.5), con condiciones iniciales no negativas, son no negativas para todo t≥0. Demostración.- Veamos que cualquier solución que entra en el cuadrante Ω, Ω := {(x, y)∈R2:x≥0, y ≥0}, permanece en él para siempre. Por la continuidad de las soluciones del sistema (1.5), cada solución debe tomar el valor cero antes que alcanzar cualquier valor negativo. Ahora bien, si tenemos y= 0 para algún instante de tiempo t=t1>0, en virtud de la segunda ecuación de (1.5), obtenemos y′(t1) = aU(x(t1))✟✟✟✯0 y(t1)−D✟✟✟✯0 y(t1) = 0, con lo que yes no decreciente en el instante de tiempo t=t1; es por esto que ynunca tomará valores negativos. De forma análoga, si x= 0 para algún instante de tiempo t=t2>0, gracias a la primera ecuación de (1.5), obtenemos x′(t2) = DI −D✟✟✟✯0 x(t2)−aU(✟✟✟✯0 x(t2))y(t2)>0, debido a que U(0) = 0,y≥0. De esta forma hemos demostrado que las soluciones no atravesarán las semirrectas positivas de los ejes de ordenadas ni abscisas; veamos ahora qué ocurre en el origen. Para ello, supongamos que existe un instante de tiempo t=t3>0tal que x(t3) = y(t3) = 0; entonces, utilizando las ecuaciones de (1.5) se llega a x′(t3) = DI −D✟✟✟✯0 x(t3)−aU(✟✟✟✯0 x(t3))✟✟✟✯0 y(t3), y′(t3) = aU(✟✟✟✯0 x(t3))✟✟✟✯0 y(t3)−D✟✟✟✯0 y(t3), de donde se deduce x′(t3)>0,y′(t3) = 0. Así, podemos concluir que las soluciones del sistema (1.5) son no negativas para todo instante de tiempo t > 0. A continuación, podemos preguntarnos el motivo por el que pueden las poblaciones de microorganismos crecer bajo un aporte limitado de nutriente, así como si el crecimiento que pueden experimentar estos microorganismos puede ser alarmante. Para responder a estas cuestiones presentamos el siguiente Teorema.
12 CAPÍTULO 1. MODELO CON UN MICROORGANISMO. Teorema 1.4.2 Cualquier solución positiva de (1.5) está acotada para condiciones iniciales no negativas que sean no nulas. Demostración.- Consideramos la función Vdada por V(t)≡V(x(t), y(t)) = x(t) + y(t). Es claro que V(0,0) = 0 yV(x(t), y(t)) ≥0, para todo instante de tiempo t≥0. Además, V(t)−→ ∞ cuando x(t), y(t)−→ ∞. Derivando Va lo largo de soluciones del sistema (1.5), obtenemos V′(t) = DI −Dx(t)−aU(x(t))y(t) + aU(x(t))y(t)−Dy(t) = =DI −D(x(t) + y(t)) =DI −DV (t). De esta manera, obtenemos V′(t) + DV (t) = DI. Para resolver esta ecuación diferencial lineal de primer orden, multiplamos por un término exponencial conveniente obteniendo eDtV′(t) + DeDtV(t) = DIeDt. Así, llegamos a que V(t) = (V(0) −I)e−Dt +I, ∀t≥0. Observar que si V(0) < I, entonces V(t)≤I,∀t≥0. Por otro lado, si V(0) > I, entonces V(t)≤V(0),∀t≥0. Así, vemos que en cualquier caso la función Vestá acotada, con lo que, por definición de V,xeytamién lo están. Fuera de la región acotada por los ejes de coordenadas positivas y la recta Dx +Dy =DI, o equivalentemente x+y=I, tenemos que V′es definida negativa, debido a que en dicha zona se verifica V=x+y > I, con lo que tenemos además que Ves una función de Lyapunov. 1.5. Puntos estacionarios o de equilibrio. El estudio de los puntos estacionarios, también llamados puntos de equilibrio, de un sistema tiene una gran importancia, que reside en que dichos puntos representan estados de energía
1.5. PUNTOS ESTACIONARIOS O DE EQUILIBRIO. 13 mínima del sistema que se esté estudiando. Los equilibrios del sistema (1.5) vienen dados por los puntos del plano xy tales que satisfacen dx dt = 0,dy dt = 0. De esta manera, los equilibrios del sistema (1.5) vienen dados por las soluciones del siguiente sistema algebraico: DI −Dx −aU(x)y= 0, aU(x)y−Dy = 0.(1.7) Una solución inmediata del sistema es (I, 0), que es un equilibrio axial. Sin embargo, si buscamos equilibrios positivos (x∗, y∗)de (1.7) deberíamos resolver el sistema DI −Dx∗−aU(x∗)y∗= 0, aU(x∗)−D= 0.(1.8) De la segunda ecuación se deduce fácilmente U(x∗) = D a, y sustituyendo en la primera ecuación de (1.8), tenemos D(I−x∗)−Dy∗= 0, de donde se obtiene I−x∗−y∗= 0, o lo que es lo mismo, x∗+y∗=I. En resumen, las condiciones que hemos obtenido son las siguientes: U(x∗) = D ayI=x∗+y∗.(1.9) Pues bien, de (1.9) podemos deducir que para obtener un equilibrio positivo, el nivel de consumo debe alcanzar, al menos una vez, el valor D/a. Esto es, debe existir un valor x∗∈(0,∞)tal que U(x∗) = D/a.
20 CAPÍTULO 1. MODELO CON UN MICROORGANISMO. 4D aβ ≥4D ae β. En conclusión, tenemos la siguiente inclusión 0,4D ae β⊂0,4D aβ, por lo que debe cumplirse 0,4D ae β. De esta forma, como o bien −V1≤V′≤ −V2, o bien −V2≤V′≤ −V1, y además hemos visto que tanto V1como V2son definidas positivas con el parámetro αadecuado, se tiene que V1>0 yV2>0, o equivalentemente −V1<0y−V2<0. Por tanto, o bien −V1≤V′≤ −V2<0, o bien −V2≤V′≤ −V1<0, es decir, V′<0, con lo que Ves definida negativa y deducimos que (x∗, y∗)es un equilibrio globalmente asintóticamente estable. Debido a que se observa, en algunos experimentos llevados a cabo en laboratorios, que el consumo sigue claramente una respuesta funcional específica dada por U(x) = x m+x, siendo m > 0la constante de saturación del medio, usaremos de aquí en adelante, salvo mención explícita que indique lo contrario, tal respuesta funcional. Haciendo uso directo de dicha función de consumo, podemos escribir el sistema (1.5) como sigue dx dt =DI −Dx −ax m+xy, dy dt =ax m+xy−Dy. (1.13) Un equilibrio positivo de (1.13) viene dado por2 U(x∗) = x∗ m+x∗=D a, 2Esto se debe al estudio hecho previamente para una función de consumo genérica. En este punto realmente lo único que se hace es particularizar lo anterior a la función de consumo que estamos tratando.
1.6. ANÁLISIS DE ESTABILIDAD. 21 es decir, x∗ m+x∗=D a, con lo que ax∗=Dm +Dx∗, es decir, ax∗−Dx∗=Dm, por tanto, (a−D)x∗=Dm, obteniendo x∗=Dm a−D. De esta forma, usando x∗+y∗=I, llegamos a Dm a−D+y∗=I, de donde y∗=I−Dm a−D. En resumen, el equilibrio que obtenemos es (x∗, y∗) = Dm a−D, I −Dm a−D. Por tanto, un conjunto de condiciones necesarias y suficientes para la existencia de equilibrio positivo para (1.13) es a(⋆) > D yI(⋆⋆) > x∗. Ahora bien, de (⋆⋆)se tiene I > mD a−D, de donde se deduce que mD I< a −D,
22 CAPÍTULO 1. MODELO CON UN MICROORGANISMO. es decir, mD I+D < a, con lo que mD +ID I< a, es decir D(m+I) I< a, y tomando inversos I D(m+I)>1 a, llegando a aI m+I (⋆⋆⋆) > D. Hemos obtenido una tercera condición, sin embargo, debemos notar que dicha condición implica (⋆), debido a que I m+I<1, con lo que D < a I m+I< a. En resumen, contábamos con dos condiciones, (⋆)y(⋆⋆); demostramos una condición equivalente a (⋆⋆), que denominamos (⋆ ⋆ ⋆)y, de ésta última, deducíamos (⋆). Por tanto, realmente la condición necesaria y suficiente para la existencia de equilibrios positivos del sistema (1.13) viene dada por D(⋆⋆⋆) < a I m+I. De esta forma, los Teoremas 1.6.1 y 1.6.2 adquieren la forma siguiente. Teorema 1.6.3 Si a≤Do bien a > D ymD a−D≥I, entonces l´ım t→∞ x(t) = Iyl´ım t→∞ y(t) = 0.
1.6. ANÁLISIS DE ESTABILIDAD. 23 Demostración.- Consideramos V(x, y) = x+yy denotamos V(t)≡V(x(t), y(t)) = x(t) + y(t). Claramente, V(0,0) = 0 yV(x, y)>0para todos x > 0,y > 0. Derivando a lo largo de soluciones del sistema tenemos V′(t) = x′(t) + y′(t) = =DI −Dx −✟✟✟✟ ✟ ax m+xy+✟✟✟✟ ✟ ax m+xy−Dy = =DI −Dx −Dy =DI −D(x+y) = DI −DV (t). Así hemos llegado a obtener la ecuación diferencial ordinaria lineal de primer orden V′(t) = −DV (t) + DI, con lo que la solución viene dada por V(t) = (V(0) −I)e−Dt +I, tal y como demostramos en el Teorema 1.6.1. Así, l´ım t→∞ V(t) = l´ım t→∞[x(t) + y(t)] = I. Como veremos que l´ım t→∞ y(t) = 0, se tendrá I= l´ım t→∞ V(t) = l´ım t→∞ x(t) +✟✟✟✟ ✟ l´ım t→∞ y(t) = l´ım t→∞ x(t), con lo que habremos concluído la demostración. Ahora procedemos distinguiendo dos casos. En primer lugar, supongamos que se verifica a≤D. En tal caso, dy dt =ax m+xy−Dy =ax m+x−Dy≤(a−D)y≤0, ya que, por hipótesis, a≤D. Ahora bien, de la última expresión se obtiene y(t)≤Ce(a−D)t,
24 CAPÍTULO 1. MODELO CON UN MICROORGANISMO. con Cconstante, debido a que tenemos una sencilla inecuación diferencial ordinaria. De esta forma, se tiene 0≤y(t)≤Ce(a−D)t−→ 0, cuando t−→ ∞, pues a−D≤0, con lo que l´ım t→∞ y(t) = 0, quedando demostrado este caso. Sin embargo, en el caso a > D yI≤mD a−Del razonamiento que debemos llevar a cabo es diferente. Sabemos que existen equilibrios si, y sólo si, se verifica la condición D < a I m+I, o equivalentemente, D < aI m+I, de donde D(m+I)< aI si, y sólo si Dm +DI < aI, es decir, Dm < aI −DI si, y sólo si Dm < I (a−D) |{z } >0 , en otras palabras, I > Dm a−D.(1.14) También sabemos que x(t) + y(t) = V(t) = (V(0) −I)e−Dt +I, con lo que l´ım t→∞[x(t) + y(t)] = I. Demostramos V′(t) = D(I−V(t)),
1.6. ANÁLISIS DE ESTABILIDAD. 25 con lo que debemos considerar dos zonas de estudio en este caso, tal y como se muestra en la siguiente figura. b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b bb x+y=I b (I, 0) x(t) y(t) x+y=I b (I, 0) x(t) y(t) u0 β0 b b Figura 1.1: Distinción de zonas. En la zona rayada, se verifica V(t) = x(t) + y(t)> I, con lo que V′<0, por tanto Ves definida negativa en tal zona. Gracias a que, además, no se verifica (1.14), no existirán equilibrios positivos, por lo que, como l´ım t→∞[x(t) + y(t)] = I, obligatoriamente las soluciones que comienzan con un dato inicial en la zona rayada, convergen al equilibrio axial (I, 0). Por otra parte, en la zona punteada partimos de datos iniciales (x0, y0)verificando x0+y0< I, de donde deducimos que x0< I; pero sabemos que I≤mD a−D, con lo que se tiene x0< I ≤mD a−D;
26 CAPÍTULO 1. MODELO CON UN MICROORGANISMO. ahora bien, y′(0) = ax0 m+x0 y0−Dy0=y0ax0 m+x0−D. Por otro lado, ax0 m+x0−D < 0si, y sólo si ax0< D(m+x0), es decir, ax0< Dm +Dx0,si, y sólo si (a−D)x0< Dm, lo que es equivalente a x0<Dm a−D, y ésto se verifica. Así, llegamos a y′(0) <0, con lo que ycomienza decreciendo; en conclusión, si el dato inicial se encuentra en la zona punteada, la solución converge al equilibrio axial (I, 0), pues en caso de crecer y atravesar la zona rayada, nos reduciríamos al caso previamente estudiado. Nótese que este argumento está basado sólidamente en la unicidad de solución del sistema. En efecto, si comienzo con un dato inicial u0en la zona punteada y la solución atraviesa la zona rayada, la solución que comenzó en u0se encuentra ahora, tras un tiempo t, en lo que puede tomarse como dato inicial para la zona rayada β0=u(t=t0; 0, u0)y entonces la solución que obtenemos tomando este dato inicial debe ser la misma que la que teníamos tomando el dato inicial u0. Por tanto, el equilibrio axial (I, 0) es un equilibrio asintóticamente estable; de esta forma, l´ım t→∞ x(t) = Iyl´ım t→∞ y(t) = 0. Nótese que las condiciónes del Teorema 1.6.3 implican la no existencia de equilibrios positivos. Teorema 1.6.4 Si a > D ymD a−D< I, entonces l´ım t→∞ x(t) = x∗=mD a−Dyl´ım t→∞ y(t) = y∗=I−x∗.
1.6. ANÁLISIS DE ESTABILIDAD. 27 Demostración.- Consideramos la función Vdada por V(t)≡V(x(t), y(t)) = 1 2(I−x−y)2+αy−y∗−y∗log y y∗. Claramente V(x∗, y∗) = 1 2(I−x∗−y∗)2+αy∗−y∗−y∗log y∗ y∗= 0 yV(x, y)>0para todos x > x∗,y > y∗, con lo que Ves una función de Lyapunov. Así, derivando Va lo largo de las soluciones de (1.13) V′= (I−x−y)(−x′−y′) + αy′−y∗1 yy′ =−(I−x−y)DI −Dx −ax m+xy+ax m+xy−Dy +αaxy m+x−Dy −y∗ax m+x+Dy∗(1.15) =−(I−x−y)DI −Dx −✟✟✟✟ ✟ ax m+xy+✟✟✟✟ ✟ ax m+xy−Dy +α(y−y∗)ax m+x−D(1.16) =−D(I−x−y)2+α(y−y∗)ax m+x−D =−D(x∗+y∗−x−y)2+α(y−y∗)ax m+x−D =−D([x−x∗] + [y−y∗])2+α(y−y∗)ax m+x−D= (⋆) Ahora bien, sabemos que x∗=mD a−D, con lo que (a−D)x∗=mD entonces ax∗−Dx∗=mD, de donde ax∗=D(m+x∗)por tanto D=ax∗ m+x∗.
28 CAPÍTULO 1. MODELO CON UN MICROORGANISMO. De esta forma, continuamos la cadena de igualdades que dejamos en (⋆)como V′=−D([x−x∗] + [y−y∗])2+α(y−y∗)axm +✘✘ ✘ axx∗−amx∗−✘✘ ✘ axx∗ (m+x)(m+x∗) =−D([x−x∗] + [y−y∗])2+αam (m+x)(m+x∗)(y−y∗)(x−x∗) Ahora bien, sabemos que D=ax∗ m+x∗entonces D x∗=a m+x∗. Continuando la cadena de igualdades V′=−D([x−x∗] + [y−y∗])2+αD x∗ m m+x(x−x∗)(y−y∗) =−D(x−x∗)2+ (y−y∗)2+ 2(x−x∗)(y−y∗)−αD x∗ m m+x(x−x∗)(y−y∗) =−D(x−x∗)2+ (y−y∗)2+2−α x∗ m m+x(x−x∗)(y−y∗). La última expresión es de la forma AX2 1+BX1Y1+CY 2 1=: H(X1, Y1), con lo que H(X1, Y1)será definida positiva si, y sólo si, A > 0yB2<4AC. De esta forma, para que Vsea definida negativa debe ocurrir 2−α x∗ m m+x2 <4si y sólo si 2−α x∗ m m+x2 <4. Pues bien, desarrollando el cuadrado tenemos 4 + α2 x∗2 m2 (m+x)2−4α x∗ m m+x<4, lo que es equivalente a 4 + α2m2−4αmx∗(m+x) x∗2(m+x)2<4si y sólo si α2m2−4αmx∗(m+x) x∗2(m+x)2<0, es decir, α2m2−4αmx∗(m+x)<0si y sólo si α2m2<4αmx∗(m+x),
1.7. EXISTENCIA DE ATRACTOR GLOBAL. 29 en otras palabras, αm < 4x∗(m+x)si y sólo si αm x∗(m+x), o lo que es lo mismo, α x∗ m m+x<4. De esta forma, V′será definida negativa y el resultado está demostrado. 1.7. Existencia de atractor global. Consideramos el sistema (1.5) dado por dx(t) dt =DI −Dx(t)−aU(x(t))y(t). dy(t) dt =aU(x(t))y(t)−Dy(t). En el resto de la sección denotaremos por u(t= 0) = u0a un dato inicial del sistema (1.5). Por otro lado, denotaremos por u(t;u0) = (x(t;t0), y(t;t0)) a una solución del sistema (1.5) y por S(t)al sistema dinámino (o semigrupo de operadores) definido por S(t)u0=S(t)(x0, y0) = (x(t), y(t)) ≡u(t;u0). Pues bien, consideramos B⊂R2un conjunto acotado. De esta forma, existe M > 0tal que B⊂B(0; M), donde B(0; M)denota la bola de centro el origen y radio Men R2. Así, para todo dato inicial u0= (x0, y0)∈B, se tiene que |u0|2≤M2. Definimos entonces w(t) := x(t) + y(t); por tanto, w′(t) = x′(t) + y′(t).
36 CAPÍTULO 2. MODELO CON DOS MICROORGANISMOS. Demostración.- (i) Consideramos la función V, dada por V(x, y1, y2) = x+y1+y2y denotamos V(t)≡V(x(t), y1(t), y2(t)) = x(t) + y1(t) + y2(t); Claramente, V(0,0,0) = 0 yV(x, y1, y2)>0para x > 0,y1>0,y2>0. Derivando Va lo largo de soluciones del sistema, tenemos V′(t) = x′(t) + y′ 1(t) + y′ 2(t) =DI −Dx −✘✘✘✘✘✘ ✘ a1 x m1+xy1−✘✘✘✘✘✘ ✘ a2 x m2+xy2+✘✘✘✘✘✘ ✘ a1 x m1+xy1−Dy1+✘✘✘✘✘✘ ✘ a2 x m2+xy2−Dy2 =DI −Dx −Dy1−Dy2 =DI −D(x+y1+y2) =DI −DV (t). Así, llegamos a V′(t) = −DV (t) + DI, con lo que V(t) = (V(0) −I)e−Dt +I, tal y como vimos en la demostración del Teorema 1.6.1. De esta forma, l´ım t→∞ V(t) = l´ım t→∞[x(t) + y1(t) + y2(t)] = I. Como posteriormente veremos que l´ım t→∞ y1(t) = 0 = l´ım t→∞ y2(t), se tiene I= l´ım t→∞ V(t) = l´ım t→∞ x(t) +✘✘✘✘✘ ✘✿0 l´ım t→∞ y1(t) +✘✘✘✘✘ ✘✿0 l´ım t→∞ y2(t) = l´ım t→∞ x(t), con lo que habremos acabado.
2.3. ANÁLISIS DE ESTABILIDAD. 37 Ahora tenemos que distinguir dos situaciones. La primera situación es la que tiene como hipótesis ai≤D,i= 1,2. En este caso, para i= 1,2, se tiene dyi dt =ai x mi+xyi−Dyi=yiaix mi+x−D≤yi(ai−D)≤0, gracias a la hipótesis y a que x mi+x≤1,i= 1,2. De esta forma yi(t)≤Ce(ai−D)t, con lo que 0≤yi(t)≤Ce(ai−D)t−→ 0, cuando t→ ∞, por lo que l´ım t→∞ yi(t) = 0, i = 1,2. Una segunda situación es aquella en que las hipótesis son ai> D yI < m´ın{λ1, λ2}. En este caso sabemos que existen equilibrios positivos si, y sólo si, se verifican ai> D yI > m´ın{λ1, λ2}(2.2) y sabemos que x(t) + y1(t) + y2(t) = V(t) = (V(0) −I)e−Dt +I, con lo que l´ım t→∞[x(t) + y1(t) + y2(t)] = I. Como ya demostramos que V′(t) = D(I−V(t)), debemos considerar dos zonas de estudio en este caso, tal y como se muestra en la siguiente figura.
38 CAPÍTULO 2. MODELO CON DOS MICROORGANISMOS. b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b bb x+y=I b (I, 0,0) x(t) (y1(t), y2(t)) x+y=I b (I, 0,0) x(t) (y1(t), y2(t)) u0 β0 b b Figura 2.1: Distinción de zonas. En la zona rayada, se verifica V(t) = x(t) + y1(t) + y2(t)> I, con lo que V′<0, por tanto Ves definida negativa en tal zona. Gracias a que, además, no se verifica (2.2), no existirán equilibrios positivos, por lo que, como l´ım t→∞[x(t) + y1(t) + y2(t)] = I, obligatoriamente las soluciones que comienzan con un dato inicial en la zona rayada, convergen al equilibrio axial (I, 0,0). Por otra parte, en la zona punteada partimos de datos iniciales (x0, y10, y20)verificando x0+y10 +y20 < I, de donde deducimos que x0< I; pero sabemos que I≤miD ai−D, i = 1,2, con lo que se tiene x0< I ≤miD ai−D, i = 1,2;
2.3. ANÁLISIS DE ESTABILIDAD. 39 ahora bien, para i= 1,2, se tiene y′ i(0) = a1 x0 mi+x0 y0i−Dy0i=y0iai x0 mi+x0−D. Por otro lado, para i= 1,2, tenemos ai x0 mi+x0−D < 0si y sólo si aix0< D(mi+x0), es decir, aix0< Dmi+Dx0,si y sólo si (ai−D)x0< Dmi, lo que es equivalente a x0<Dmi ai−D, y ésto se verifica. Así, llegamos a y′ i(0) <0,i= 1,2, con lo que yi,i= 1,2, comienza decreciendo, luego la suma y1+y2también decrece; en conclusión, si el dato inicial se encuentra en la zona punteada, la solución converge al equilibrio axial (I, 0,0), pues en caso de crecer y atravesar la zona rayada, nos reduciríamos al caso previamente estudiado. Conviene notar que este argumento está basado sólidamente en la unicidad de solución del sistema. En efecto, si comienzo con un dato inicial u0en la zona punteada y la solución atraviesa la zona rayada, la solución que comenzó en u0se encuentra ahora, tras un tiempo t, en lo que puede tomarse como dato inicial para la zona rayada β0=u(t=t0; 0, u0)y entonces la solución que obtenemos tomando este dato inicial debe ser la misma que la que teníamos tomando el dato inicial u0. Por tanto, el equilibrio axial (I, 0,0) es un equilibrio asintóticamente estable y de esta forma, l´ım t→∞ x(t) = Iyl´ım t→∞ y1(t) = 0 = l´ım t→∞ y2(t). (ii) En este caso, para i= 1,2, se tiene ai x mi+xyi−Dyi=ai x mi+xyi−Dmi+x mi+xyi= (ai−D)x mi+xyi−miD mi+xyi=
40 CAPÍTULO 2. MODELO CON DOS MICROORGANISMOS. = (ai−D)x mi+xyi−(ai−D)λi yi mi+x= (ai−D)x mi+xyi−miDyi mi+x= = (ai−D)x mi+xyi−λi yi mi+x= (ai−D)x−λi mi+xyi. Así, (2.1) queda dx dt =DI −Dx −a1 x m1+xy1−a2 x m2+xy2, dy1 dt = (a1−D)x−λ1 m1+xy1,(2.3) dy2 dt = (a2−D)x−λ2 m2+xy2. Consideramos la función de Lyapunov V, dada por V(t)≡V(x(t), y1(t), y2(t)) = x−λ1−λ1log x λ1 +c1y1−y∗ 1−y∗ 1log y1 y∗ 1+c2y2, con ci=ai ai−D,i= 1,2, constantes positivas. Por tanto, V′(t) = x′−λ1 1 ✚ ✚ λ1 x ✚ ✚ λ1 x′+c1y′ 1−y∗ 1 1 y1 y′ 1+c2y′ 2 =x−λ1 xDI −Dx −a1 x m1+xy1−a2 x m2+xy2 +c11−y∗ 1 y1(a1−D)x−λ1 m1+xy1+c2y2(a2−D)x−λ2 m2+x =x−λ1 xDI −Dx −a1 x m1+xy1−a2 x m2+xy2 +c1 y1−y∗ 1 ✚ ✚ y1 (a1−D)x−λ1 m1+x✚ ✚ y1+c2y2(a2−D)x−λ2 m2+x =x−λ1 x(DI −Dx)−a1 x−λ1 m1+xy1−a2 x−λ2 m2+xy2 +a1y1 x−λ1 m1+x−a1y∗ 1 x−λ1 m1+x+a2 x−λ2 m2+xy2 =x−λ1 x(DI −Dx)−a1y∗ 1 x−λ1 m1+x+a2 λ1−λ2 m2+xy2 = (x−λ1)D(I−x) x−a1 y∗ 1 m1+x+a2 λ1−λ2 m2+xy2.
2.3. ANÁLISIS DE ESTABILIDAD. 41 Consideramos ahora el equilibrio (x∗=λ1, y∗ 1,0). Por (2.3), entonces dx dt =DI −Dx −a1 x m1+xy1−a2 x m2+xy2, y evaluando en el equilibrio, dx∗ dt =DI −Dx∗−a1 x∗ m1+x∗y∗ 1−✘✘✘✘✘✘✘ ✘✿0 a2 x∗ m2+x∗·0, es decir, dλ1 dt =DI −Dλ1−a1 λ1 m1+λ1 y∗ 1⇔0 = D(I−λ1)−a1 λ1 m1+λ1 y∗ 1, de donde a1 λ1 m1+λ1 y∗ 1=D(I−λ1), o equivalentemente, y∗ 1=D(I−λ1)m1+λ1 a1λ1 . Con esto, el término D(I−x) x−a1 y∗ 1 m1+x=D(I−x) x−D(I−λ1)(m1+λ1) λ1(m1+x) =D(I−x)λ1(m1+x)−Dx(I−λ1)(m1+λ1) λ1x(m1+x) =D λ1x(m1+x){(I−x)λ1(m1+x)−x(I−λ1)(m1+λ1)} =D λ1x(m1+x)λ1Im1+✟✟ ✟ λ1Ix −✘✘✘ ✘ λ1xm1−λ1xx −xIm1 −✟✟ ✟ xIλ1+✘✘✘ ✘ xλ1m1+xλ1λ1 =D λ1x(m1+x){−m1I(x−λ1)−λ1x(x−λ1)} =−D(x−λ1) λ1m(m1+x){m1I−λ1x}.
42 CAPÍTULO 2. MODELO CON DOS MICROORGANISMOS. Así, V′(t) = −D(x−λ1)2 λ1x(m1+x){m1I+λ1x}+a2 λ1−λ2 m2+xy2≤0, pues λ1< λ2. De esta forma, V′es semidefinida negativa. A continuación vamos a demostrar que M:= {(x∗=λ1, y∗ 1,0)} es el mayor conjunto invariante de E:= {(x∗=λ1, y1,0), y1≥0}. y1(t) y2(t) x(t) + λ1E b ×M y∗ 1 Figura 2.2: Conjuntos MyE. Consideramos el conjunto Edefinido previamente. Ahora bien, dicho conjunto es realmente E={(x∗=λ1, y1,0), y1≥0}.
2.3. ANÁLISIS DE ESTABILIDAD. 43 Por otro lado x′E=DI −Dλ1−a1 λ1 m1+λ1 y10. Por tanto, x′E= 0 si y sólo si DI −Dλ1−a1 λ1 m1+λ1 y10 = 0, o lo que es lo mismo D(I−λ1) = a1 λ1 m1+λ1 y10 si y sólo si D(m1+λ1)(I−λ1) a1λ1 , y esto ocurre si y sólo si y10 =y∗ 1. Esto quiere decir (obsérvese Figura 2.2) que si empiezo con un dato inicial distinto de y∗ 1, me salgo del conjunto E(este caso es el que reflejan las trayectorias punteadas en la Figura 2.2), con lo que el conjunto Manteriormente definido es el mayor subconjunto invariante de E. En conclusión, (x∗, y∗ 1,0) es asintóticamente estable, gracias a una consecuencia del Teorema de La Salle (véase ??). Notar que la demostración de la segunda parte es exactamente la misma, intercambiando los papeles de λ1yλ2. (iii) Extendiendo la función de Lyapunov usada en el Teorema 1.6.4, consideramos V(t)≡V(x(t), y1(t), y2(t)) = 1 2(I−x−y1−y2)2+α 2 X i=1 yi−y∗ i−y∗ ilog yi y∗ i. Además, V(x∗, y∗ 1, y∗ 2) = 1 2(I−x∗−y∗ 1−y∗ 2)2+α 2 X i=1 y∗ i−y∗ i−y∗ ilog y∗ i y∗ i =1 2(I−x∗−y∗ 1−y∗ 2)2=1 2(y∗ 1+y∗ 2−y∗ 1−y∗ 2)2= 0. Además, V(x, y1, y2)>0para x > x∗,y1> y∗ 1,y2> y∗ 2. Así, derivando Va lo largo de soluciones de (2.1) y procediendo como en el Teorema 1.6.4, tenemos
44 CAPÍTULO 2. MODELO CON DOS MICROORGANISMOS. V′(t) = (I−x−y1−y2)(−x′−y′ 1−y′ 2) + α 2 X i=1 y′ i−y∗ i y′ i y∗ i yi y∗ i = (I−x−y1−y2)(−x′−y′ 1−y′ 2) + α 2 X i=1 y′ i−y∗ i y′ i yi = (I−x−y1−y2)−DI +Dx +a1 x m1+xy1+a2 x m2+xy2 +(D−a1)x−λ1 m1+xy1+ (D−a2)x−λ2 m2+xy2 +α 2 X i=1 (ai−D)x−λi mi+xyi−y∗ i(ai−D)x−λi mi+x✓✓✓ yi 1 yi = (I−x−y1−y2)(−DI +Dx +a1 x m1+xy1+a2 x m2+xy2+Dx−λ1 m1+xy1 −a1 x−λ1 m1+x+Dx−λ2 m2+xy2−a2 x−λ2 m2+xy2 +α 2 X i=1 (ai−D)x−λi mi+xyi−y∗ i(ai−D)x−λi mi+x = (I−x−y1−y2)−DI +Dx +✘✘✘✘✘✘ ✘ a1 x m1+xy1+✘✘✘✘✘✘ ✘ a2 x m2+xy2+Dx−λ1 m1+xy1 +Dx−λ2 m2+xy2+a1 λ1 m1+xy1+a2 λ2 m2+xy2−✘✘✘✘✘✘ ✘ a1 x m1+xy1−✘✘✘✘✘✘ ✘ a2 x m2+xy2 +α 2 X i=1 (yi−y∗ i)(ai−D)x−λi mi+x= (∗). Ahora bien, sabemos que λi=miD ai−D,i= 1,2, por tanto Dx m1+xy1+Dx m2+xy2−Dλ1 m1+xy1−Dλ2 m2+xy2+a1 λ1 m1+xy1+a2 λ2 m2+xy2 =Dx m1+xy1+Dx m2+xy2+ (a1−D)λ1 m1+xy1+ (a2−D)λ2 m2+xy2 = 2 X i=1 Dx mi+xyi+ (ai−D)miD (mi+x)(ai−D)yi
2.3. ANÁLISIS DE ESTABILIDAD. 45 = 2 X i=1 Dx mi+xyi+Dmi mi+xyi =D 2 X i=1 yi✚✚✚✚ x+mi x+mi =D[y1+y2]. Entonces V′= (I−x−y1−y2)(−DI +Dx +D[y1+y2]) + α 2 X i=1 (yi−y∗ i)(ai−D)x−λi mi+x. Nos centramos en la expresión α 2 X i=1 (yi−y∗ i)(ai−D)x−λi mi+x =α 2 X i=1 (yi−y∗ i)(ai−D)x mi+x−(yi−y∗ i)(ai−D)λi mi+x =α 2 X i=1 ((yi−y∗ i)(ai−D)x mi+x−(yi−y∗ i)(ai−D) Dmi ai−D mi+x) =α 2 X i=1 (yi−y∗ i)(ai−D)x mi+x−(yi−y∗ i)✘✘✘✘ ✘ (ai−D)Dmi ✘✘✘✘ ✘ (ai−D)(mi+x) =α 2 X i=1 (yi−y∗ i)aix−Dx −Dmi mi+x =α 2 X i=1 (yi−y∗ i)aix mi+x−Dx mi+x−Dmi mi+x ai−D ai−D =α 2 X i=1 (yi−y∗ i)aix mi +x−Dx mi +x−(ai−D)x∗ mi+x =α 2 X i=1 (yi−y∗ i)(x−x∗)ai−D mi+x= (⋆).
52 CAPÍTULO 3. MODELO CON CRECIMIENTO EN PARED. cuando el ratio de flujo es bajo, la biomasa muerta no se elimina del sistema inmediatamente y está sujeta a descomposición bacteriana, lo que a su vez conduce a la regeneración del nutriente. No se espera que el 100 % de la biomasa sea reciclada, sino sólamente una fracción que denotaremos por b∈(0,1). Cuando IyDson ambas constantes y no hay retardos en el sistema, como es el caso, el modelo que se presenta a continuación describe la dinámica de quimiostatos con crecimiento en pared. Observemos que sólo y1(t)contribuye al material reciclado de la biomasa muerta en el medio. Además, como los microorganismos adheridos a la pared no se limpian del sistema, el término −Dy2(t)no se incluye en la ecuación que representa el crecimiento de y2(t). Todos los parámetros son los mismos que aparecen en (1.5), es decir, el modelo con un microorganismo que estudiamos en el Capítulo 1, a diferencia de que 0< c ≤areemplaza a acomo el ratio de crecimiento de las especies consumidoras1. De esta forma, obtenemos el siguiente modelo dx dt =DI −Dx −ax m+xy1−ax m+xy2+bνy1, dy1 dt =−(ν+D)y1+cx m+xy1−r1y1+r2y2,(3.1) dy2 dt =−νy2+cx m+xy2+r1y1−r2y2. 3.2. Puntos estacionarios o de equilibrio. En este caso, los equilibrios del sistema (3.1) vienen dados por las soluciones de DI −Dx −ax m+xy1−ax m+xy2+bνy1= 0, −(ν+D)y1+cx m+xy1−r1y1+r2y2= 0,(3.2) −νy2+cx m+xy2+r1y1−r2y2= 0. Claramente (I, 0,0) es solución inmediata del sistema. Sin embargo, si intentamos seguir los pasos que llevamos a cabo a la hora de estudiar los puntos de equilibrio de (3.1), el modelo con dos microorganismos del Capítulo 2, y buscamos equilibrios de la forma (x∗, y∗ 1,0) o(x∗,0, y∗ 2), nos damos cuenta de que, en este caso, tales equilibrios no son posibles. Veamos una explicación de este hecho. Supongamos que, por ejemplo, y2= 0, entonces de la tercera ecuación del sistema (3.2) tenemos r1y1= 0, siendo r1>0, con lo que se deduce que y1= 0 y, por tanto, x=I; en resumen, 1Notemos que hemos unificado las constantes respecto al modelo estudiado en el Capítulo 2, modelo con dos microorganismos, para hacer más sencillo el estudio del mismo de aquí en adelante.
3.2. PUNTOS ESTACIONARIOS O DE EQUILIBRIO. 53 obtenemos el equilibrio (I, 0,0). Si, por el contrario, tenemos y1= 0, de la segunda ecuación del sistema (3.2) tenemos r2y2= 0, siendo r2>0, con lo que se deduce que y2= 0 y, por tanto, x=I; en resument, volvemos a obtener el equilibrio (I, 0,0). Por tanto, además del equilibrio trivial (I, 0,0) lo que debemos buscar ahora son equilibrios positivos. Sin embargo, antes de continuar, vamos a realizar un cambio de variables en el sistema (3.1), que viene dado por α(t) = y1(t) y1(t) + y2(t), z(t) = y1(t) + y2(t). Así, el sistema resultante es el siguiente dx dt =DI −Dx −ax m+xz+bναz, dz dt =−νz −Dαz +cx m+xz, (3.3) dα dt =−Dα(1 −α)−r2α+r2(1 −α). En efecto, es trivial la primera ecuación que se obtiene tras el cambio. En cuanto a la segunda z′=y′ 1+y′ 2 =−νy1−Dy1+cx m+xy1−✘✘ ✘ r1y1+✘✘ ✘ r2y2−νy2+cx m+xy2+✘✘ ✘ r1y1−✘✘ ✘ r2y2 =−ν(y1+y2)−Dy1+cx m+x(y1+y2) =−νz −Dy1+cx m+xz. Ahora bien, por definición de α, se tiene y1=α(y1+y2) = αz, por tanto, se deduce que z′=−νz +cx m+xz−Dαz.
54 CAPÍTULO 3. MODELO CON CRECIMIENTO EN PARED. Finalmente, para la tercera ecuación α′=y1 y1+y2′ =y′ 1(y1+y2)−y1(y′ 1+y′ 2) (y1+y2)2 =✟✟ ✟ y′ 1y1+y′ 1y2−✟✟ ✟ y1y′ 1−y1y′ 2 (y1+y2)2 =y′ 1y2−y1y′ 2 (y1+y2)2. Por otro lado y′ 1y2=−(ν+D)y1y2 + cx m+xy1y2−r1y1y2+r2y2 2, −y1y′ 2=νy2y1−cx m+xy2y1−r1y2 1+r2y2y1. De aquí, como α′=y′ 1y2−y1y′ 2 (y1+y2)2, se tiene y′ 1y2−y1y′ 2= (y1+y2)2α′, con lo que (y1+y2)2α′=−Dy1y2−r1(y1y2+y2 1) + r2(y2 2+y2y1). En conclusión α′=−Dy1y2 (y1+y2)2−r1 y1y2+y2 1 (y1+y2)2+r2 y2 2+y2y1 (y1+y2)2. Ahora bien α(1 −α) = y1 y1+y21−y1 y1+y2=y1 y1+y2 y2 y1+y2 =y1y2 y1+y2 ; y1y2+y2 1 (y1+y2)2=y1✘✘✘✘ ✘ (y2+y1) (y1+y2)✁ 2=y1 y1+y2 =α; y2 2+y2y1 (y1+y2)2=y2✘✘✘✘ ✘ (y2+y1) (y1+y2)✁ 2=y2 y1+y2 = 1 −α.
3.2. PUNTOS ESTACIONARIOS O DE EQUILIBRIO. 55 Con todos los cálculos realizados previamente llegamos a la ecuación diferencial para α α′=−Dα(1 −α)−r1α+r2(1 −α), tal y como nosotros queríamos, quedando así justificadas cada una de las ecuaciones que forman el sistema (3.3). Para cualesquiera soluciones positivas y1,y2de (3.3) se tiene 0< α(t)<1, por la propia definición de α, para todo t. Además α′α=0 = [−Dα(1 −α)−r1α+r2(1 −α)] α=0 =r2>0, α′α=1 = [−Dα(1 −α)−r1α+r2(1 −α)] α=1 =−r1<0. En consecuencia, el intervalo (0,1) es una región positivamente invariante para α(t). Por otro lado, la ecuación que verifica αes independiente de xy de zy sus puntos de equilibrio son las soluciones de la siguiente ecuación de segundo grado −Dα(1 −α)−r1α+r2(1 −α) = 0, es decir, Dα2+ (−r2−D−r1)α+r2= 0, y que resulta ser α∗ ±=D+r1+r2±p(D+r1+r2)2−4Dr2 2D, donde α∗ +denota la elección de signo positivo en el numerador y α∗ −la elección de signo negativo. Antes de continuar veamos que α∗ −es la única de las dos raíces que pertenece al intervalo (0,1) y que α∗ +es mayor que 1. En efecto, (D+r1+r2)2−4Dr2=D2+r2 1+r2 2−2r1r2+ 2Dr1+ 2Dr2−4Dr2 =D2+r2 1+r2 2+ 2r1r2+ 2Dr1−2Dr2> D2+r2 1+r2 2+ 2r1r2−2Dr1−2Dr2 = (D−(r1+r2))2= (D−r1−r2)2, con lo que p(D+r1+r2)2−4Dr2>p(D−r1−r2)2=p(r1+r2−D)2,
56 CAPÍTULO 3. MODELO CON CRECIMIENTO EN PARED. de donde se deduce p(D+r1+r2)2−4Dr2>p(r1+r2−D)2=r1+r2−D, obteniendo −p(D+r1+r2)2−4Dr2<−(r1+r2−D) = D−r1−r2, es decir D+r1+r2−p(D+r1+r2)2−4Dr2< D +r1+r2+D−r1−r2= 2D. En consecuencia α∗ −=D+r1+r2−p(D+r1+r2)2−4Dr2 2D<2D 2D= 1. Ahora demostraremos que α∗ +es mayor que 1. En efecto, p(D+r1+r2)2−4Dr2<p(D+r1+r2)2=D+r1+r2, de donde se deduce −p(D+r1+r2)2−4Dr2>−D−r1−r2, llegando a D+r1+r2−p(D+r1+r2)2−4Dr2> D +r1+r2−D−r1−r2= 0. En consecuencia, α∗ −>0 2D= 0. Que α∗+>1se deduce inmediatamente pues D+r1+r2+p(D+r1+r2)2−4Dr2 2D>1si y sólo si p(D+r1+r2)2−4Dr2> D−r1−r2,
3.2. PUNTOS ESTACIONARIOS O DE EQUILIBRIO. 57 lo que es equivalente a (D+r1+r2)2−4Dr2>(D−r1−r2)2 y esto es cierto. De aquí en adelante notaremos α∗a la elección negativa, es decir, denotaremos así a α∗ −. A continuación veamos que α(t)−→ α∗cuando t−→ +∞y, por tanto, α∗es asintóticamente estable. En efecto, α(t)verifica la ecuación diferencial ordinaria α′(t) = −Dα(1 −α)−r1α+r2(1 −α) =−Dα +Dα2−r1α+r2−r2α =Dα2−(D+r1+r2)α+r2. A continuación, utilizando el cambio α(t) = α∗+ 1/z(t), ya que α∗es solución, obtenemos z′(t) = (D+r1+r2−2Dα∗)z−D, que tiene como solución z=e(D+r1+r2−2Dα∗)tC+D D+r1+r2−2Dα∗. Deshaciendo el cambio inicial obtenemos α(t) = α∗+D+r1+r2−2Dα∗ (D+r1+r2−2Dα∗)Ce(D+r1+r2−2Dα∗)t+D y tomando α(t) = α∗se deduce que γ γCeγt∗+D= 0,donde γ=D+r1+r2−2Dα∗, por tanto, como γ6= 0, es claro que t∗= +∞, donde t∗denota el instante de tiempo para el que se alcanza α∗. Sustituyendo αpor α∗en el sistema (3.3) tras eliminar la tercera ecuación del mismo tenemos dx dt =DI −Dx −ax m+xz+bνα∗z, dz dt =−νz −Dα∗z+cx m+xz. (3.4)
58 CAPÍTULO 3. MODELO CON CRECIMIENTO EN PARED. Los equilibrios del sistema (3.4) vienen dados por la solución de DI −Dx∗−ax∗ m+x∗z∗+bνα∗z∗= 0, −νz∗−Dα∗z∗+cx∗ m+x∗z∗= 0.(3.5) De la segunda ecuación de (3.5) obtenemos −ν−Dα∗+cx∗ m+x∗z∗= 0. Ahora bien, por definición de α(t)yz(t),z∗= 0 no puede ser solución del sistema (3.5), con lo que −ν−Dα∗+cx∗ m+x∗= 0, o, en otras palabras, x∗ m+x∗=ν+Dα∗ c De aquí podemos obtener la expresión explicítica de x∗ x∗ m+x∗=ν+Dα∗ c⇔cx∗=mν +mDα∗+νx∗+Dα∗x∗⇔ ⇔(Dα∗+ν−c)x∗+mν +mDα∗= 0 ⇔(c−ν−Dα∗)x∗=m(ν+Dα∗), es decir x∗=m(ν+Dα∗) c−ν−Dα∗.(3.6) Por otro lado DI −Dx∗+z∗−ax∗ m+x∗+bνα∗= 0, o, lo que es lo mismo, z∗=DI −Dx∗ ax∗ m+x∗−bνα∗.(3.7) De esta forma hemos obtenido el equilibrio del sistema (3.4) y observándo las expresiones obtenidas es claro que, para que existan equilibrios positivos, se deben verificar ν+Dα∗< c y(I−x∗)ax∗ m+x∗−bνα∗>0,
3.3. ANÁLISIS DE ESTABILIDAD. 59 o, equivalentemente, ν+Dα∗< c y(I−x∗)ax∗ m+x∗−bνα∗>0.(3.8) 3.3. Análisis de estabilidad. Estudiaremos en esta sección la estabilidad de los puntos de equilibrio obtenidos en la sección previa. Recogemos toda la información en el siguiente Teorema. Teorema 3.3.1 Se verifican los siguientes resultados: (i) Si ν+Dα∗> c, entonces el equilibrio axial (I, 0) del sistema (3.4) es globalmente asintóticamente estable, es decir, cualquier solución de (3.4) satisface l´ım t→∞ x(t) = I, l´ım t→∞(y1(t) + y2(t)) = l´ım t→∞ z(t) = 0. (ii) Si ν+Dα∗< c,D−c+az∗>0yx∗< I, entonces el equilibrio positivo (x∗, z∗)del sistema (3.4) es globalmente asintóticamente estable, es decir, cualquier solución de (3.4) satisface l´ım t→∞ x(t) = x∗,l´ım t→∞(y1(t) + y2(t)) = l´ım t→∞ z(t) = z∗. Demostración.- (i) Consideremos la segunda ecuación del sistema (3.4) dz dt =−νz −Dα∗z+cx m+xz=−ν+Dα∗−cx m+xz. Sabemos que x m+x≤1,entonces −cx m+x≥ −c. De esta forma, se tiene ν+Dα∗−cx m+x≥ν−Dα∗−c.
60 CAPÍTULO 3. MODELO CON CRECIMIENTO EN PARED. Así, como ν+Dα∗> c, se tiene ν+Dα∗−c > 0, con lo que ν+Dα∗−cx m+x≥ν−Dα∗−c > 0. Se obtiene por tanto dz dt =−ν+Dα∗−cx m+x<0, lo que implica que zes decreciente; además, z(t) = y1(t) + y2(t)≥0, por lo que l´ım t→∞(y1(t) + y2(t)) = l´ım t→∞ z(t) = 0. Por otro lado, usando la primera ecuación del sistema (3.4) se tiene dx dt =DI −Dx −ax m+xz+bνα∗z. Ahora bien, dado ε > 0, existe T > 0tal que ∀t≥Tse verifica que −ax(t) m+x(t)+bνα∗z(t)≤ε. Entonces consideramos dx dt ≤DI −Dx +ε, ∀t≥T, lo que nos da x(t)≤x(T)e−Dt +I+ε D y de dx dt ≥DI −Dx +ε se deduce x(t)≤x(T)e−Dt +I−ε D.
3.3. ANÁLISIS DE ESTABILIDAD. 61 Entonces l´ım t−→∞ x(t)∈hI−ε D, I +ε Di y como ε > 0es arbitrario, tenemos l´ım t−→∞ x(t) = I. (ii) Para llevar a cabo el estudio de la estabilidad del equilibrio positivo (x∗, z∗)vamos a reescribir el sistema (3.4), usando (3.5), como sigue dx dt =−Dx +Dx∗−ax m+x(z−z∗)−az∗x m+x−x∗ m+x∗+bνα∗(z−z∗), dz dt =cz x m+x−x∗ m+x∗.(3.9) En efecto, partimos de (3.5), es decir DI −Dx∗−ax∗ m+x∗z∗+bνα∗z∗= 0, −νz∗−Dα∗z∗+cx∗ m+x∗z∗= 0. De la primera ecuación se tiene DI =Dx∗+ax∗ m+x∗z∗−bνα∗z∗, con lo que, sustituyendo en la primera de las ecuaciones, se obtiene dx dt =Dx∗+ax∗ m+x∗z∗−bνα∗z∗−Dx −ax m+xz+bνα∗z. Sumando y restando ax m+xz∗
68 CAPÍTULO 3. MODELO CON CRECIMIENTO EN PARED. Pues bien, definimos ahora bkcomo bk:= m´ın{D, k}. De esta forma w′≤DI −bkw. Además, sabemos que w(0) = x(0) + y1(0) + y2(0) = x0+y10 +y20. Resolviendo el problema de Cauchy dw dt ≤ −bkw +DI, w(0) = w0=x0+y10 +y20, obtenemos una cota superior para w(t); veámoslo. Consideramos la EDO dw dt ≤ −bkw +DI y la escribiremos equivalentemente como dw dt +bkw ≤DI. Multiplicamos por eb kt y obtenemos eb kt dw dt +bkeb ktw≤DIeb kt, es decir, d dt heb ktwi≤DIeb kt. Integrando llegamos a que eb ktw(t)−w(t0)≤Zt t0 DIeb ksds =DI bkeb kt −eb kt0≤DI bkeb kt −1, pues t0≥0, con lo que −eb kt0≤ −1. También podemos escribir lo anterior, equivalentemente, como
3.4. EXISTENCIA DE ATRACTOR GLOBAL. 69 w(t)≤e−b ktw(t0) + DI bk1−e−b kt =e−b ktw(t0)−DI bke−b kt +DI bk =w(t0)−DI bke−b kt +DI bk =x0+y10 +y20 −DI bke−b kt +DI bk. En consecuencia, w(t)≤Ce−b kt +DI bk, C =x0+y10 +y20 −DI bk. De esta manera, usando una vez más la conocida desigualdad de Young y que x(t), y1(t), y2(t)≥ 0, llegamos a |S(t)u0|2=x2(t) + y2 1(t) + y2 2(t) ≤ |x(t) + y1(t) + y2(t)|2 =|w(t)|2 ≤Ce−b kt +DI bk2 =(x0+y10 +y20)e−b kt +DI bk−DI bke−b kt2 =(x0+y10 +y20)e−b kt +DI bk1−e−b kt2 ≤2(x0+y10 +y20)2e−2b kt + 2 DI bk21−e−b kt |{z } ≤1 2. Por otro lado, también sabemos que (x0+y10 +y20)2≤3|(x0, y10, y20)|2, tal y como demostramos en la Sección 2.4.
70 CAPÍTULO 3. MODELO CON CRECIMIENTO EN PARED. Por tanto, |S(t)u0|2≤6|(x0, y10, y20)|2e−2b kt + 2 DI bk2 ≤6M2e−2b kt + 2 DI bk2 ≤1 + 2 DI bk2 =: ρ2, pues l´ım t→∞ 6M2e−2b kt = 0. De esta forma, existe TB≥0tal que para todo t≥TB, se tiene 6M2e−2b kt ≤1; en efecto, 6M2e−2b kt ≤1⇔6M2≤e2b kt ⇔log(6M2)≤2bkt ⇔1 2bklog(6M2)≤t, por tanto, definimos TB:= 1 2bklog(6M2). En resumen, hemos demostrado que existe un compacto, B(0, ρ), de tal forma que existe TB≥0 tal que para todo t≥TB, se tiene |S(t)u0| ≤ ρ; por tanto, existe un compacto absorbente, con lo que existe un atractor global que denotaremos por A.
Parte II Modelo no autónomo sin retardo. 71
Capítulo 4 Ratio de consumo de nutriente variable. 4.1. Introducción. En esta parte del trabajo vamos a comenzar estudiando lo que ocurre cuando algunas de las constantes de los modelos que estudiábamos en la primera parte se vuelven dependientes de la variable tiempo. Comenzamos este capítulo estudiando el caso en el que la concentración de nutriente es constante pero el ratio de consumo del mismo es variable. Específicamente vamos a asumir que Dvaría de forma continua en tiempo, por ejemplo, de forma periódica o aleatoria, en un intervalo acotado positivo D(t)∈(dm, dM)para todo t∈R. En el resto del capítulo trataremos en primer lugar el caso del modelo sin pared para, posteriormente, introducir los nuevos cambios al modelo en el que permitimos crecimiento en pared. 4.2. Modelo sin pared. En esta sección estudiamos el caso del modelo sin pared. Cuando Ies una constante positiva y Dvaría en tiempo, tenemos el sistema dx(t) dt =D(t)(I−x(t)) −ax(t) m+x(t)y(t), dy(t) dt =−D(t)y(t) + ax(t) m+x(t).(4.1) Lema 4.2.1 Para cualquier tiempo inicial t0∈Ry cualesquiera condiciones iniciales x0, y0≥ 0, todas las soluciones del sistema (4.1) son no negativas y acotadas para todo t≥t0. 73
74 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. Demostración.- Los coeficientes son continuamente diferenciables para x, y ≥0. En particular, el término no lineal axy m+x=ay 1−m m+x es no negativo y acotado por la función lineal ay en el cuadrante positivo. En efecto, aym son constantes postivias y x, y ≥0, con lo que dicho término no lineal es no negativo; además, sabemos que x m+x≤1, con lo que axy m+x≤ay en el primer cuadrante, ya que x, y ≥0. Esto asegura la existencia y unicidad de soluciones para cualquier tiempo para el que permanezcan en el primer cuadrante. En efecto, en virtud de la segunda ecuación del sistema (4.1) y del hecho de que x(t) m+x(t)≤1, tenemos dy(t) dt =−D(t)y(t) + ax(t) m+x(t)y(t) ≤ −D(t)y(t) + ay(t) = (a−D(t))y(t). De esta forma, resolviendo la inecuación diferencial que nos resulta tenemos y(t)−y(t0)≤Zt t0 (a−D(s))y(s)ds ≤(a−dm)Zt t0 y(s)ds, por tanto y(t)≤y(t0) + (a−dm)Zt t0 y(s)ds. Aplicando el Lema de Gronwall, se tiene 0≤y(t)≤y(t0)e(a−dm)(t−t0), por tanto y(t)está definida para cualquier instante de tiempo, con lo que efectivamente tenemos existencia de solución. La unicidad se tiene gracias a los resultados de la teoría básica de ecuaciones diferenciales ordinarias.
4.2. MODELO SIN PARED. 75 Por la continuidad de las soluciones, con condición inicial x(t0) = x0>0,x(t)tiene que tomar el valor 0antes de tomar cualquier valor negativo. Pero dx dt x=0 =D(t)I > 0, con lo que x(t)no puede tomar valores negativos. Con la condición inicial y(t0) = y0>0, existe t1> t0tal que y(t)>0en [t0, t1]. Además, y(t) = y0eRt t0(−D(s)+ax(s) m+x(s))ds para t∈[t0, t1]. Debido a la unicidad de solución, esta expresión se cumple para todo t≥t0, así y(t)es no negativa. Sumando las ecuaciones del sistema (4.1) tenemos d(x(t) + y(t)) dt =−D(t)(x(t) + y(t)−I) de donde se deduce inmediatamente que cuando x(t) + y(t)≥I, entonces I≤x(t) + y(t)≤ x0+y0. De forma análoga, cuando x(t) + y(t)≤I, tendremos 0≤x(t) + y(t)≤I. Por tanto, 0≤x(t) + y(t)≤m´ax{x0+y0, I}, con lo que x(t)yy(t)están acotadas. Ahora pasamos a estudiar el comportamiento a largo plazo de las soluciones del sistema (4.1). Más específicamente, estableceremos condiciones bajo las cuales el atractor es una solución completa o un único punto. Notemos que (I, 0) es el único equilibrio del sistema para todos los valores de los parámetros. Otras soluciones atrayentes no serán puntos de equilibrio. Teorema 4.2.1 Supongamos que D:R−→ [dm, dM], donde 0< dm< dM<∞, es continua. Entonces el sistema (4.1) tiene un atractor pullback A={A(t) : t∈R}en el cuadrante positivo R2 +:= {(x, y)∈R2:x≥0, y ≥0}.Además, (i) Si a < dm, el equilibrio axial (I, 0) es asintóticamente estable en el cuadrante positivo y el atractor pullback Atiene una única componente que es el subconjunto A(t) = {(I, 0)} para todo t∈R; (ii) si lo que sucede es que a > 1 + m IdM
76 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. entonces el atractor pullback Atambién contiene puntos estrictamente dentro del cuadrante positivo, además del punto {(I, 0)}; (iii) cuando, además de la hipótesis de (ii) se tiene dm< a < dm(mdm+dMI)2 (mdm+dMI)2−mId2 m el atractor pullback Aestá formado por el punto de equilibrio axial {(I, 0)}y una solución completa ξ∗que está uniformemente acotada fuera de los ejes, así como soluciones completas heteroclínicas entre ellas, es decir, sus componentes son A(t) = {(x, y)∈R2 +:x+y=I;ξ∗(t)≤x≤I} para todo t∈R. Demostración.- Definiendo w(t) = x(t) + y(t)y sumando las ecuaciones del sistema (4.1) llegamos a dw(t) dt =dx(t) dt +dy(t) dt =D(t)I−D(t)x(t)−✘✘✘✘✘✘✘ ✘ ax(t) m+x(t)y(t)−D(t)y(t) +✘✘✘✘✘✘✘ ✘ ax(t) m+x(t)y(t) =D(t)I−D(t)[x(t) + y(t)] =D(t)[I−w(t)]. Esta ecuación diferencial ordinaria posee un equilibrio w∗que viene dado por la solución de 0 = D(t)[I−w∗], es decir, w∗=I, incluso cuando D(t)no es constante. Se puede demostrar que es atrayente, tanto en el sentido pullback como en el sentido forward, veámoslo. Sea w0:= w(t0) = x(t0) + y(t0), donde x(·),y(·)denotan la solución de (4.1), con datos iniciales (t0; (x0(t0), y0(t0))). Entonces, resolviedo la EDO dw(t) dt =D(t)(I−w(t))
4.2. MODELO SIN PARED. 77 obtenemos w(t) = w0e−Rt t0D(s)ds +Ie−Rt t0D(s)ds Zt t0 D(s)eRs t0D(r)drds =w0e−Rt t0D(s)ds +Ie−Rt t0D(s)ds Zt t0 d ds heRs t0D(r)drids =w0e−Rt t0D(s)ds +Ie−Rt t0D(s)dseRt t0D(r)dr =w0e−Rt t0D(s)ds +I, que converge a Ien las dos siguientes situaciones: (i) cuando t0−→ −∞, estando tfijo, y (ii) t−→ ∞, estando t0fijo, ya que 0≤e−Rt t0D(s)ds ≤e−dm(t−t0)−→ 0 en ambos casos. En virtud del Lema 4.2.1, se tiene que para todo ε > 0, el conjunto compacto no vacío Bε:= {(x, y)∈R2 +:x+y≤I+ε} es positivamente invariante y absorbente en R2 +. Por tanto, el sistema dinámico no autónomo en R2 +generado por el sistema de ecuaciones diferenciales ordinarias (4.1) tiene un atractor pullback A={A(t) : t∈R} formado por subconjuntos compactos no vacíos de R2 +.
84 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. De esta forma, por una parte, suponiendo aI > (I+m)dM, tenemos que dmI−dMx(t)−aI < dmI−dMx(t)−(I+m)dM =−(dM−dm)I−dMx(t)−mdM<0. Y, por otra parte, dMI−dmx(t) = dM(I−x(t)) + (dM−dm)x(t)>0. Entonces, para cualesquiera dos soluciones x1(t)yx2(t)de (4.2), ∆(t) := x1(t)−x2(t) satisface d∆(t) dt =dx1(t) dt −dx2(t) dt = (I−x1(t)) D(t)−ax1(t) m+x1(t)−(I−x2(t)) D(t)−ax2(t) m+x2(t) =✟✟✟ ✟ ID(t)−x1(t)D(t)−aI x1(t) m+x1(t)+ax2 1(t) m+x1(t)−✟✟✟ ✟ ID(t) + x2(t)D(t) +aI x2(t) m+x2(t)−ax2 2(t) m+x2(t) =−D(t)(x1(t)−x2(t)) −aI x1(t)(m+x2(t)−x2(t)(m+x1(t)) (m+x1(t))(m+x2(t)+ +ax2 1(t)(m+x2(t)) −x2 2(t)(m+x1(t)) (m+x1(t))(m+x2(t)) =−D(t)(x1(t)−x2(t)) −aI mx1(t) +✘✘✘✘✘ ✘ x1(t)x2(t)−mx2(t)−✘✘✘✘✘ ✘ x1(t)x2(t) (m+x1(t))(m+x2(t)) +amx2 1(t) + x2 1(t)x2(t)−mx2 2(t)−x1(t)x2 2(t) (m+x1(t))(m+x2(t)) =−D(t)∆(t)−aIm (m+x1(t))(m+x2(t))∆(t) + am(x1(t) + x2(t)) + x1(t)x2(t) (m+x1(t))(m+x2(t)) ∆(t), ya que (m(x1(t) + x2(t)) + x1(t)x2(t))∆(t) = (mx1(t) + mx2(t) + x1(t)x2(t))(x1(t)−x2(t)) =mx2 1(t)−✘✘✘✘✘✘ ✘ mx1(t)x2(t) +✘✘✘✘✘✘ ✘ mx1(t)x2(t)−mx2 2(t) + x2 1(t)x2(t)−x1(t)x2 2(t) =mx2 1(t)−mx2 2(t) + x2 1(t)x2(t)−x1(t)x2 2(t).
4.2. MODELO SIN PARED. 85 Por (4.4), obtenemos, para i= 1,2, xi(t)≤dMI dm , entonces m+xi(t)≤m+dMI dm , por tanto, (m+x1(t))(m+x2(t)) ≤m+dMI dm2 y, por ende, −amI (m+x1(t))(m+x2(t)) ≤ − amI m+dMI dm2. Además, m(x1(t) + x2(t)) + x1(t)x2(t) = mx1(t) + mx2(t) + x1(t)x2(t) < m2+mx1(t) + mx2(t) + x1(t)x2(t) = (m+x1(t))(m+x2(t)), con lo que m(x1(t) + x2(t)) + x1(t)x2(t) (m+x1(t))(m+x2(t)) <1. Así, d∆(t) dt =−D(t)∆(t)−amI (m+x1(t))(m+x2(t))∆(t) +am(x1(t) + x2(t)) + x1(t)x2(t) (m+x1(t))(m+x2(t)) ∆(t) <−dm∆(t)−amI m+dMI dm2∆(t) + a∆(t),
86 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. con lo que ∆(t)−→ 0cuando t−→ ∞, siempre y cuando dm+amI m+dMI dm2> a, es decir, cuando a < dm(mdm+dMI)2 (mdm+dMI)2−mId2 m . En efecto, a < dm+amI m+dMI dm2=dm+amI m2+d2 MI2 d2 m+ 2mdMI dm =dm+amI d2 mm2+d2 MI2+2mdMdmI d2 m =dm+amId2 m d2 mm2+d2 MI2+ 2mdMdmI; entonces a−amId2 m d2 mm2+d2 MI2+ 2mdMdmI< dm, es decir, a1−mIdm2 d2 mm2+d2 MI2+ 2mdMdmI< dm, con lo que ad2 mm2+d2 MI2+ 2mdMdmI−mId2 m d2 mm2+d2 MI2+ 2mdMdmI< dm, por tanto, se concluye que a < dm(d2 mm2+d2 MI2+ 2mdMdmI) d2 mm2+d2 MI2+ 2mdMdmI−mId2 m ,
4.3. MODELO CON PARED. 87 o lo que es lo mismo, a < dm(mdm+dMI)2 (mdm+dMI)2−mId2 m . Esto se cumple si a < dmcomo en el caso (1). Sin embargo, puede también cumplirse si aes ligeramente mayor que dm. En este caso el límite pullback para condiciones iniciales estrictamente positivas del sistema escalar (4.2) es uniformemente estrictamente contractivo [38],[39] en (0, I)y existe una solución completa ξ∗(t)∈(0, I), que es también asintóticamente estable en el sentido forward usual. El correspondiente atractor pullback A1de este sistema en [0, I]incluye el equilibrio Iy los conjuntos A1(t) = [ξ∗(t), I], para cada t∈R, es decir, incluye las trayectorias heteroclínicas que unen los equilibrios ξ∗(t) yI. Para el sistema 2−dimensional (4.1) el atractor pullback Aestá formado por los conjuntos A(t) = {(x, y) : x+y=I;ξ∗(t)≤x≤I} en R2 +, para t∈R. 4.3. Modelo con pared. Estudiaremos en esta sección lo que ocurre cuando hacemos variar el ratio de consumo de nutriente en el modelo con pared. Cuando Ies constante y Dvaría en tiempo, el modelo que tenemos es el siguiente dx(t) dt =D(t)(I−x(t)) −ax(t) m+x(t)(y1(t) + y2(t)) + bνy1(t), dy1(t) dt =−(ν+D(t))y1(t) + cx(t) m+x1(t)y1(t)−r1y1(t) + r2y2(t),(4.5) dy2(t) dt =−νy2(t) + cx(t) m+x(t)y2(t) + r1y1(t)−r2y2(t), donde recordamos que aes el ratio máximo de crecimiento específico, crepresenta el ratio de crecimiento de las especies consumidoras, con lo que a≥c;mes la constante de saturación del medio de consumo; r1, r2representan los ratios con los que las especies pasan de estar en el medio a estar en la pared, o viceversa; νdenota el ratio de muerte colectiva de y; y bdescribe la fracción de biomasa muerta que es reciclada.
88 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. Como las variables x, y1, y2representan concentraciones, es claro que suponemos condiciones iniciales no negativas: x(t0) = x0;y1(t0) = y10;y2(t0) = y20. Lema 4.3.1 Supongamos que (x0, y10, y20)∈R3 +:= {(x, y1, y2)∈R3:x≥0, y1≥0, y2≥ 0}. Entonces todas las soluciones del sistema (4.5) correspondientes a datos iniciales en R3 + son (i) no negativas para todo t > t0; (ii) uniformemente acotadas en R3 +. Además, el sistema dinámico no autónomo sobre R3 +generado por el sistema de ecuaciones diferenciales ordinarias (4.5) tiene un atractror pullback A={A(t) : t∈R}en R3 +. Demostración.- (i) Por continuidad, cada solución tiene que tomar el valor 0antes de tomar cualquier valor negativo. Con x= 0,y1≥0,y2≥0, la ecuación diferencial ordinaria para x(t)se reduce a dx(t) dt =D(t)I+bνy1(t)>0, con lo que x(t)es estrictamente creciente en x= 0. Con y1= 0,x≥0,y2≥0, la EDO reducida para y1(t)es dy1(t) dt =r2y2≥0, con lo que y1(t)es no decreciente en y1= 0. De forma similar, con y2= 0,x≥0,y1≥0, la EDO reducida para y2(t)es dy2(t) dt =r1y1(t)≥0, con lo que y2(t)también es no decreciente en y2= 0. Además, (x(t), y1(t), y2(t)) ∈R3 + para cualquier instante de tiempo t. (ii) Definimos, para X(t) = (x(t), y1(t), y2(t)) ∈R3 +, kX(t)k1:= x(t) + y1(t) + y2(t).
4.3. MODELO CON PARED. 89 Entonces kX(t)k1≤S(t)≤a ckX(t)k1, donde S(t) = x(t) + a c(y1(t) + y2(t)). En efecto, sabemos que 0< c ≤a, con lo que a/c ≥1. De esta forma, kX(t)k1=x(t) + y1(t) + y2(t)≤x(t) + a c(y1(t) + y2(t)) =S(t) = x(t) + a c(y1(t) + y2(t)) ≤a c(x(t) + y1(t) + y2(t)) = a ckX(t)k1. Derivando S(t)a lo largo de soluciones del sistema (4.5) tenemos dS(t) dt =dx(t) dt +a c dy1(t) dt +a c dy2(t) dt =D(t)I−D(t)x(t)−✘✘✘✘✘✘✘✘✘✘✘✘ ✘ ax(t) m+x(t)(y1(t) + y2(t)) + bνy1(t)−a c(ν+D(t))y1(t) +✘✘✘✘✘✘✘ ✘ ax(t) m+x(t)y1(t)−✟✟✟✟ ✟ r1a cy1(t) +✟✟✟✟ ✟ r2a cy2(t)−aν cy2(t) +✘✘✘✘✘✘✘ ✘ ax(t) m+x(t)y2(t) +✟✟✟✟ ✟ a cr1y1(t)−✟✟✟✟ ✟ a cr2y2(t) =D(t)[I−x(t)] −ha c(ν+D(t)) −bνiy1(t)−a cνy2(t) ≤dMI−dmx(t)−ha c(ν+dm)−bνiy1(t)−a cνy2(t).(4.6) Notamos que 0< c ≤ay0< b < 1, entonces a c(ν+dm)−bν ≥ν+dm−bν > ν +dm−ν=dm. Sea µ:= m´ın{dm, ν}, entonces dS(t) dt ≤dMI−dmx(t)−ha c(ν+dm)−bνiy1(t)−a cνy2(t) < dMI−dmx(t)−dmy1(t)−νy2(t) ≤dMI−µ{x(t) + y1(t) + y2(t)} =dMI−µS(t).
90 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. De esta forma llegamos a dS(t) dt < dMI−µS(t).(4.7) Tratando esta inecuación exactamente como hicimos en la demostración de (iii)del Teorema 4.2.1, llegamos a d dt eµtS(t)< dMIeµt. Integrando entre t0ytobtenemos eµtS(t)−S(t0)<dMI µeµt −eµt0 <dMI µeµt −1, dado que t0≥0. Si continuamos realizando cuentas tal como en la demostración de (iii) del Teorema 4.2.1, llegamos a S(t)< S(t0)e−µt +dMI µ(1 −e−µt) |{z } ≥0 . Así, si S(t0)<dMI µ, entonces S(t)<dMI µe−µt +dMI µ(1 −e−µt) < ✚✚✚✚ ✚ dMI µe−µt +dMI µ−✚✚✚✚ ✚ dMI µe−µt =dMI µ,∀t≥t0. Por otra parte, si S(t0)≥dMI µ, entonces S(t)< S(t0)e−µt +S(t0)(1 −e−µt) =✘✘✘✘✘ S(t0)e−µt +S(t0)−✘✘✘✘✘ S(t0)e−µt =S(t0),∀t≥t0.
4.3. MODELO CON PARED. 91 En resumen, si S(t0)<dMI µ, entonces S(t)<dMI µpara todo t≥t0. Por otro lado, si S(t0)≥dMI µ, entonces S(t)será no creciente para todo t≥t0, con lo que S(t)< S(t0). Entonces implica que kX(t)k1está acotada superiormente, es decir, kX(t)k1≤m´ax dMI µ, x(t0) + a c(y1(t0) + y2(t0)), para todo t≥t0. Se sigue entonces que para todo ε > 0, el conjunto compacto no vacío Bε:= (x, y1, y2)∈R3 +:x+a c(y1+y2)≤dMI µ+ε es positivamente invariante y absorbente en R3 +. El sistema dinámico no autónomo en R3 +generado por el sistema de ecuaciones diferenciales ordinarias (4.5) tiene por tanto un atractor pullback A={A(t) : t∈R}, que consiste en subconjuntos compactos no vacíos de R3 +tal que están contenidos en Bε. Para obtener más información sobre la estructura interna del atractor pullback del sistema no autónomo generado por el sistema de ecuaciones diferenciales ordinarias (4.5), hacemos el cambio de variables α(t) = y1(t) y1(t) + y2(t), z(t) = y1(t) + y2(t),(4.8) tal y como hicimos en el estudio que realizamos para la primera parte del trabajo. De esta forma, y análogamente a lo que hicimos allí, obtenemos el sistema dx(t) dt =D(t)(I−x(t)) −ax(t) m+x(t)z(t) + bνα(t)z(t), dz(t) dt =−νz(t)−D(t)α(t)z(t) + cx(t) m+x(t)z(t),(4.9) dα(t) dt =−D(t)α(t)(1 −α(t)) −r1α(t) + r2(1 −α(t)). Observemos que el punto de equilibrio (I, 0,0) del sistema (4.5) no cuenta para el sistema (4.9), ya que αno está definido para él. Por otra parte, (I, 0) es un equilibrio para las dos primeras ecuaciones del sistema (4.9).
92 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. 4.3.1. Dinámica global de α(t). Observamos que la dinámica de α(t) = α(t, t0, α0)está desacoplada de x(t)yz(t), y satisface la ecuación de Ricatti que aparece como tercera ecuación del sistema (4.9). Para cualesquiera y1,y2positivos, tenemos 0< α(t)<1para todo t. Recordemos que α′α=0 =r2>0yα′α=1 =−r1<0, con lo que el intervalo (0,1) es positivamente invariante, tal y como estudiamos con anterioridad en el trabajo. Ésta es la región biológicamente relevante. Cuando Des constante, ya estudiamos que existe un único punto crítico α∗∈(0,1), dado por α∗:= D+r1+r2−p(D+r1+r2)2−4Dr2 2D(4.10) Queremos ahora tratar el caso en que Dvaría en tiempo, aleatoria o periódicamente, en un intervalo positivo y acotado D(t)∈[dm, dM]para todo t∈R. En este caso necesitamos hablar de atractores pullback aleatorios o determinísticos Aα={Aα(t) : t∈R}en el intervalo (0,1). Tal atractor existe, ya que el intervalo unidad es positivamente invariante, con lo que sus componentes están dadas por Aα(t) = \ t0<t α(t, t0,[0,1]),∀t∈R. Estas componentes son de la forma Aα= [α∗ l(t), α∗ u(t)], donde α∗ l(t)yα∗ u(t)son soluciones completas y acotadas para la ecuación de Ricatti. El resto de soluciones completas y acotadas se encuentran entre ellas. Podemos usar inecuaciones diferenciales para obtener cotas de estas soluciones completas. Para ello, lo primero que hacemos es reescribir la ecuación de Ricatti del sistema (4.9) de la forma α′(t) = D(t)(α2(t)−α(t)) −(r1+r2)α(t) + r2,(4.11) tal y como justificamos en estudios de capítulos previos. Como α(t)<1yD(t)>0, entonces α′(t)≤ −(r1+r2)α(t) + r2. Así, α(t)≤β(t), con α(t0) = β(t0), donde β′(t) = −(r1+r2)β(t) + r2.
4.3. MODELO CON PARED. 93 Esta EDO tiene un punto estacionario asintóticamente estable, dado por β∗=r2 r1+r2 , con lo que las soluciones completas de la ecuación de Ricatti del sistema (4.9) se encuentran (menos una cantidad infinitesimal) por debajo, es decir, α∗ u(t)≤β∗para t∈R. Esto nos proporciona una cota superior. Por otro lado, tenemos α′(t) = D(t)α2(t)−(D(t) + r1+r2)α(t) + r2 ≥ −(dM+r1+r2)α(t) + r2. Así, α(t)≥γ(t), con α(t0) = γ(t0), donde γ′(t) = −(dM+r1+r2)γ(t) + r2. Esta EDO tiene un punto crítico asintóticamente estable, dado por γ∗:= r2 r1+r2+dM . En este caso obtenemos una cota inferior α∗ l(t)≥γ∗para todo t∈R. En resumen, A(t) = [α∗ l(t), α∗ u(t)] ⊂[γ∗, β∗]. Para investigar el caso en que el atractor pullback se compone de una única solución completa, necesitamos encontrar condiciones bajo las cuales α∗ l(t)≡α∗ u(t), t ∈R. Supongamos que no son iguales y consideremos su diferencia ∆α(t) := α∗ u(t)−α∗ l(t). Entonces ∆′ α(t) = D(t)(α∗2 u(t)−α∗ u(t)) −(r1+r2)α∗ u(t) + r2 −D(t)(α∗2 l(t)−α∗ l(t)) + (r1+r2)α∗ l(t)−r2 =D(t)α∗2 u(t)−D(t)α∗ u(t)−(r1+r2)α∗ u(t) −D(t)α∗2 l(t) + D(t)α∗ l(t) + (r1+r2)α∗ l(t) =D(t)(α∗2 u(t)−α∗2 l(t)) −(D(t) + r1+r2)(α∗ u(t)−α∗ l(t)) =D(t)(α∗ u(t) + α∗ l(t))(α∗ u(t)−α∗ l(t)) −(D(t) + r1+r2)(α∗ u(t)−α∗ l(t)) =D(t)(α∗ u(t) + α∗ l(t))∆α(t)−(D(t) + r1+r2)∆α(t) ≤2dMα∗ u(t)∆α(t)−(dm+r1+r2)∆α(t) ≤2dMβ∗∆α(t)−(dm+r1+r2)∆α(t) =2dMr2 r1+r2 ∆α(t)−(dm+r1+r2)∆α(t) =2dMr2 r1+r2−dm−r1−r2∆α(t).
100 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE. En efecto, q1=dm a−c+ν+dM−bνβ∗, con lo que cq1I m+q1I= cIdm a−c+ν+dM−bνβ∗ m+dmI a−c+dM−bνβ∗ = cIdm ✭✭✭✭✭✭✭ ✭ a−c+ν+dM−bνβ∗ m(a−c+dM−bνβ∗)+dMI ✭✭✭✭✭✭✭ ✭ a−c+ν+dM−bνβ∗ =cIdm m(a−c+dM−bνβ∗) + dMI. Así, para ε > 0suficientemente pequeño dz(t) dt z=ε =−ν−D(t)α∗(t) + cx(t) m+x(t)ε >−ν−dMβ∗+c(q1I−ε) m+q2I−ε>0.(4.15) Las desigualdades (4.14), (4.15), junto con d(x(t) + z(t)) dt x+z=q1I >0yd(x(t) + z(t)) dt x+z=q2I <0,(4.16) aseguran la invariancia positiva del conjunto compacto Tεy la existencia de un atractor pullback A={A(t) : t∈R} en Tε. Vamos a demostrar, para concluir la demostración de este Teorema, las desigualdades (4.16). Por un lado
4.3. MODELO CON PARED. 101 du(t) dt u=q1I > dmI−(a−c+ν+dM−bνβ∗)q1I(4.17) =I[dm−(a−c+ν+dM−bνβ∗)q1] ="dm−✭✭✭✭✭✭✭✭✭✭✭✭ ✭ (a−c+ν+dM−bνβ∗)dm ✭✭✭✭✭✭✭✭✭✭✭ ✭ a−c+ν+dM−bνβ∗# =I[dm−dm] = 0. De esta forma, d(x(t) + z(t)) dt x+z=q1I =du(t) dt u=q1I >0. Análogamente tenemos la segunda desigualdad. du(t) dt u=q2I < dMI−dmγ∗q2I=I(dM−dmγ∗q2) =IdM−✟✟ ✟ dmγ∗dM ✟✟ ✟ dmγ∗=I(dM−dM) = 0, con lo que d(x(t) + z(t)) dt x+z=q2I =du(t) dt u=q2I <0. Desafortunadamente, en esta ocasión no somos capaces de obtener la existencia de una solución completa estable que atraiga a todas las soluciones completas positivas como se hizo en el caso del modelo sin crecimiento en pared.
102 CAPÍTULO 4. RATIO DE CONSUMO DE NUTRIENTE VARIABLE.
Capítulo 5 Ratio de concentración de nutriente variable. 5.1. Introducción. Estudiamos en este capítulo el caso en el que el ratio de consumo de nutriente es constante, pero el ratio de concentración del mismo es variable. Aquí vamos a suponer que el ratio de concentración de nutriente Ipuede variar de forma continua con el paso del tiempo, por lo que lo denotaremos por I(t). De forma similar a como hicimos en el capítulo anterior, suponemos que Iestá acotada con valores positivos, en particular, I(t)∈[im, iM]para todo t∈R, donde 0< im< iM<∞. En el resto del capítulo trataremos en primer lugar el caso del modelo sin pared para, posteriormente, introducir los nuevos cambios al modelo en el que permitimos crecimiento en pared. 5.2. Modelo sin pared. Consideramos en esta sección el caso del modelo sin pared, donde el sistema de ecuaciones diferenciales ordinarias que estudiaremos es el siguiente dx(t) dt =D(I(t)−x(t)) −ax(t) m+x(t)y(t), dy(t) dt =−Dy(t) + ax(t) m+x(t)y(t).(5.1) Sea w(t) := x(t) + y(t). Entonces 103
104 CAPÍTULO 5. RATIO DE CONCENTRACIÓN DE NUTRIENTE VARIABLE. dw(t) dt =D(I(t)−w).(5.2) Esta ecuación no tiene equilibrios cuando I(t)no es constante, pero tiene un equilibrio no autónomo no trivial1que es tanto atrayente pullback como forward: w(t) = w(t0)e−D(t−t0)+De−D(t−t0)Zt t0 I(s)eD(s−t0)ds =w(t0)e−D(t−t0)+De−Dt Zt t0 I(s)eDsds, que converge a w∗(t) = De−Dt Zt −∞ I(s)eDsds cuando t0−→ −∞ y cuando t−→ ∞, es decir, l´ım t→∞ |w(t)−w∗(t)|= l´ım t0→−∞ |w(t)−w∗(t)|= 0. Notar que w∗(t)∈[im, iM]para todo t∈Rdebido a las cotas para I. Lema 5.2.1 Para cualquier tiempo inicial t0∈Ry cualesquiera condiciones iniciales x0, y0≥ 0, todas las soluciones del sistema (5.1) son no negativas y acotadas para todo t≥t0. Demostración.- Los coeficientes son continuamente diferenciables para x, y ≥0. En particular, el término no lineal axy m+x=ay 1−m m+x es no negativo y acotado por la función lineal ay en el cuadrante positivo, tal y como vimos en la demostración del Lema 4.2.1. Esto segura la existencia y unicidad de soluciones para cualquier tiempo para el que permanezcan en el primer cuadrante. Por la continuidad de las soluciones, con condición inicial x(t0) = x0>0,x(t)tiene que tomar el valor 0antes de tomar cualquier valor negativo. Pero dx(t) dt x=0 =DI(t)>0, con lo que x(t)no puede tomar valores negativos. Con la condición inicial y(t0) = y0>0, existe t1> t0tal que y(t)>0en [t0, t1]. 1de hecho es una solución especial que se construye usando una técnica pullback.
5.2. MODELO SIN PARED. 105 Además, y(t) = y0e−D(t−t0)+Rt t0ax(s) m+x(s)ds para t∈[t0, t1]. Debido a la unicidad de solución, esta expresión se cumple para todo t≥t0, así y(t)es no negativa. Teorema 5.2.1 El sistema dinámico no autónomo generado por el sistema de ecuaciones diferenciales ordinarias (5.1) tiene un atractor pullback A={A(t) : t∈R}en R2 +. Además, (i) cuando D > a, la solución completa (x∗(t), y∗(t)) = (w∗(t),0) es asintóticamente estable en R2 +y el atractor pullback tiene una única componente A(t) = {(w∗(t),0)}para todo t∈R; (ii) cuando, además, aim> D(m+iM), el atractor pullback tiene componentes no triviales que incluyen (w∗(t),0) y puntos estrictamente positivos; (iii) cuando tenemos la hipótesis de (ii)y, además, D < a ya(m2+m(2iM−im) + i2 M)< D(m+iM)2, el atractor pullback contiene una solución completa no trivial que atrae a todas las otras soluciones completas positivas. Demostración2.- En virtud del Lema 5.2.1 y del hecho de que w∗(t)∈[im, iM], el conjunto compacto no vacío B:= {(x, y)∈R2 +:im≤x+y≤iM} es positivamente invariante y absorbente en R2 +para la EDO (5.2). 2Antes de nada, comentaremos que en la demostración de este Teorema se omitirán algunos detalles por analogía a demostraciones del capítulo anterior.
106 CAPÍTULO 5. RATIO DE CONCENTRACIÓN DE NUTRIENTE VARIABLE. x y x+y=im b b im b b x+y=iM iM Figura 5.1: Representación gráfica del conjunto B. El sistema dinámico no autónomo en R2 +generado por el sistema de ecuaciones diferenciales ordinarias (5.1) tiene así un atractor pullback A={A(t) : t∈R}formado por subconjuntos no vacíos de B. Entonces (w∗(t),0) ∈A(t)para todo t∈R, ya que el atractor pullback contiene a todas las soluciones completas acotadas. Para probar (i)notamos que la primera ecuación del sistema (5.1) se puede acotar superiormente como sigue dy(t) dt =−D−ax(t) m+x(t)y(t)≤ −(D−a)y(t), con lo que se tiene de forma inmediata que y(t)−→ 0cuando t−→ ∞, para D > a. En cuanto a (ii), del signo positivo de la derivada de la segunda ecuación del sistema (5.1), tenemos que x(t)es creciente en la zona x= 0 del conjunto absorbente B. La zona y= 0 es invariante, pero para y=ε << imyim≤x≤iM, la primera ecuación del sistema (5.1) da dy(t) dt =ax(t) m+x(t)−Dy(t)≥aim m+iM−Dy(t)>0, cuando aim> D(m+iM). Esto significa que el interior positivo del conjunto absorbente también contiene puntos del atractor pullback. Finalmente, para ver (iii)consideramos la primera EDO del sistema (5.1) restringida a la variedad estable x(t) + y(t) = w∗(t), sobre la cual adquiere la forma dx(t) dt =D(I(t)−x(t)) −ax(t) m+x(t)(w∗(t)−x(t)).(5.3)
5.2. MODELO SIN PARED. 107 Para cualesquiera dos soluciones x1(t),x2(t)de (5.3), definimos ∆x(t) := x1(t)−x2(t). Entonces ∆xsatisface d∆x(t) dt =dx1(t) dt −dx2(t) dt =✟✟✟ ✟ DI(t)−Dx1(t)−ax1(t) m+x1(t)w∗(t) + ax2 1(t) m+x1(t)−✟✟✟ ✟ DI(t) +Dx2(t) + ax2(t) m+x2(t)w∗(t)−ax2 2(t) m+x2(t) =−Dx1(t) + Dx2(t)−ax1(t)(m+x2(t)) −x2(t)(m+x1(t)) (m+x1(t))(m+x2(t)) w∗(t) +ax2 1(t)(m+x2(t)) −x2 2(t)(m+x1(t)) (m+x1(t))(m+x2(t)) =−D∆x(t)−amw∗(t) (m+x1(t))(m+x2(t))∆x(t) + am(x1(t) + x2(t)) + x1(t)x2(t) (m+x1(t))(m+x2(t)) ∆x(t), donde la última igualdad se obtiene por razonamientos idénticos a los que se realizaron en la demostración de la parte (iii)del Teorema 4.2.1. Como 0≤x(t)≤w∗(t)≤iMyw∗(t)≥im, tenemos d∆x(t) dt =−D∆x(t)−aamw∗(t) (m+x1(t))(m+x2(t))∆x(t) + am(x1(t) + x2(t)) + x1(t)x2(t) (m+x1(t))(m+x2(t)) ∆x(t) <−D∆x(t)−amim (m+iM)2∆x(t) + a∆x(t), análogamente a la demostración de la parte (iii)del Teorema 4.2.1. Así, ∆x(t)−→ 0cuando t−→ ∞, siempre y cuando D+amim (m+iM)2> a, o, equivalentemente a(m2+m(2iM−im) + i2 M)< D(m+iM)2, ya que a < D +amim (m+iM)2,
108 CAPÍTULO 5. RATIO DE CONCENTRACIÓN DE NUTRIENTE VARIABLE. con lo que a1−mim (m+iM)2< D, es decir, am2+i2 M+ 2miM−mim (m+iM)2< D, o lo que es lo mismo, a(m2+m(2iM−im) + i2 M)< D(m+iM)2. Esto también se cumple si a < D, en cuyo caso estamos en la situación (i)del Teorema. Sin embargo, puede también cumplirse si aes ligeramente mayor que (m2+m(2iM−im) + i2 M)<(m+im)2, en cuyo caso las estimaciones anteriores con x1(t),x2(t)igual a w∗(t), el sistema es contractivo uniformemente estrictamente [38], [39] en el cuadrante positivo y entonces tiene una única solución completa como su atractor pullback en el cuadrante positivo. 5.3. Modelo con pared. Estudiaremos en esta sección lo que ocurre cuando hacemos variar el ratio de concentración de nutriente en el modelo con pared. Cuando Des constante e Ivaría en tiempo, el modelo que tenemos es el siguiente dx(t) dt =D(I(t)−x(t)) −ax(t) m+x(t)(y1(t) + y2(t)) + bνy1(t), dy1(t) dt =−(ν+D)y1(t) + cx(t) m+x1(t)y1(t)−r1y1(t) + r2y2(t),(5.4) dy2(t) dt =−νy2(t) + cx(t) m+x(t)y2(t) + r1y1(t)−r2y2(t), Lema 5.3.1 Supongamos que (x0, y10, y20)∈R3 +. Entonces, todas las soluciones del sistema (5.4) con datos iniciales (x(t0), y1(t0), y2(t0)) = (x0, y10, y20)son (i) no negativas para todo t≥t0; (ii) uniformemente acotadas en R3 +.
5.3. MODELO CON PARED. 109 Además, el sistema dinámico no autónomo en R3 +generado por el sistema de ecuaciones diferenciales ordinarias (5.4) tiene un atractor pullback A={A(t) : t∈R}en R3 +. Demostración.- (i) Por continuidad, cada solución tiene que tomar el valor 0antes de tomar cualquier valor negativo. Con x= 0,y1≥0,y2≥0, la EDO para x(t)se reduce a dx(t) dt =DI(t) + bνy1>0, con lo que x(t)es estrictamente creciente en x= 0. Con y1= 0,x≥0,y2≥0, la EDO reducida para y1(t)es dy1(t) dt =r2y2≥0, por tanto, y1(t)es no decreciente en y1= 0. De forma similar, con y2= 0,x≥0,y1≥0, la EDO reducida para y2(t)es dy2(t) dt =r1y1≥0, por lo que y2(t)es no decreciente en y2= 0. Además, (x(t), y1(t), y2(t)) ∈R3 +para cualquier t∈R. (ii) Definimos, para (x(t), y1(t), y2(t)) ∈R3 +, kX(t)k1:= x(t) + y1(t) + y2(t). Entonces, kX(t)k1≤S(t)≤a ckX(t)k1, donde S(t) = x(t) + a c(y1(t) + y2(t)), análogamente a como se hacía en la demostración de la parte (ii)del Lema 4.3.1.
Capítulo 6 Modelo sin crecimiento en pared. 6.1. Introducción. Otra posibilidad a la hora de estudiar los modelos de quimiostatos es considerar un retardo en tales modelos de tal forma que quede explicado el crecimiento oscilatorio de los microorganismos. Los retardos son naturales en los sistemas biológicos, debido a que permiten la coexistencia de poblaciones en competición como una solución periódica no forzada. La presencia de términos con retardo en modelos de quimiostatos está completamente justificada, ya que el comportamiento de los sistemas dinámicos no sólo depende en general del presente, sino también de su historia. A veces sólo un corto período de tiempo pasado es el que proporciona información relevante al sistema (es el caso de retardos constantes o variables) y otras veces es toda la historia la que hay que tener en cuenta (retardos distribuídos). Debemos notar que la presencia de retardos no constantes hace que estemos trabajando con modelos de quimiostatos no autónomos. Los primeros modelos de quimiostatos con retardo fueron estudiados por Caperon en [35] Más tarde, Bush y Cook propusieron un modelo de quimiostato con retardo en el ratio de crecimiento intrínseco de los microorganismos, pero sin retardo en la ecuación del nutriente [36]. Podemos introducir entonces un retardo al modelo más simple de los que hemos estudiado dx(t) dt =DI −Dx −aU(x(t))y(t), dy(t) dt =−Dy(t) + cU(x(t))y(t). 117
118 CAPÍTULO 6. MODELO SIN CRECIMIENTO EN PARED. resultando el siguiente modelo dx(t) dt =DI −Dx −aU(x(t))y(t), dy(t) dt =−Dy(t) + cU(x(t−ρ(t)))y(t),(6.1) siendo ρ(·) : R+−→ [0, h] (h > 0) una función continua. En este caso consideramos un retardo variable, en lugar de uno constante, aunque también podría tenerse en cuenta un retardo distribuído. La elección del retardo variable es para mostrar el comportamiento de estos modelos a la hora de introducir un retardo, sin entran en terrenos que pueden ser pantanosos en el caso de retardos distribuídos y que pudieran hacer perder la visión de lo que se pretende mostrar. 6.2. Algunas consideraciones previas. Para resolver el sistema (6.1), necesitamos fijar los valores de la solución no sólo en un tiempo inicial t0∈R, sino en el intervalo [t0−h, t0], es decir, necesitamos encontrar las soluciones al sistema (6.1) para t≥t0tal que x(t) = φ1(t−t0), y(t) = φ2(t−t0),para t∈[t0−h, t0],(6.2) donde (φ1, φ2)∈Ch:= C([−h, 0]; R2). Observemos que el sistema (6.1) se puede escribir como la siguiente ecuación diferencial funcional no autónoma dz(t) dt =f(t, zt),(6.3) donde f(·,·) : R×Ch→R2está dada por f(t, φ) = DI −Dφ1(0) −aU(φ1(0))φ2(0) −Dφ2(0) + cU(φ1(−ρ(t))φ2(0) ,para φ:= φ1 φ2∈Ch, yzt∈Chdenota la función zt(θ) = z(t+θ),para θ∈[−h, 0], para cualquier función continua z: [−h, +∞)→R2. Para resolver el problema de valores iniciales para (6.3), una posibilidad es considerar una condición inicial φ∈Chen el tiempo inicial t0∈R. Esto nos proporcionará, bajo condiciones adecuadas, la única solución z(·;t0, φ)de (6.3) tal que z0(·;t0, φ) = φ, lo que significa que
6.3. SIGNO Y ACOTACIÓN DE SOLUCIONES. 119 zt0(θ;t0, φ) := z(t0+θ;t0, φ) = φ(θ)para θ∈[−h, 0]. De esta forma podemos construir un sistema dinámico no autónomo (por ejemplo, un proceso) en el espacio de fases Chen el siguiente sentido: para cualquier t≥t0, definimos el operador Z(t, t0) : Ch→Chcomo Z(t, t0)φ=zt(·;t0, φ),para φ∈Ch.(6.4) 6.3. Signo y acotación de soluciones. Al igual que en los capítulos previos, la región biológicamente interesante en nuestro modelo es aquella donde las soluciones sean positivas y finitas. En el próximo teorema demostraremos que las soluciones del sistema (6.1) son no negativas y acotadas. Lema 6.3.1 Dada cualquier condición inicial φ≥0en [t0−h, t0], las soluciones de (6.1) son no negativas. Demostración.- Lo primero que vamos a demostrar es que si una solución comienza en el cuadrante R2 +={(x, y) : x≥0, y ≥0}, permanece allí para siempre. Por continuidad, cada solución debe tomar el valor 0antes de alcanzar cualquier valor negativo. De hecho, como dy(t) dt y=0 = [−Dy(t) + cU(x(t−ρ(t)))y(t)] y=0 = 0, entonces yno decrece en y= 0, con lo que nunca tomará valores negativos. Por otro lado, dx(t) dt x=0 = [DI −Dx(t)−aU(x(t))y(t)] x=0 =DI −a✟✟✟✯0 U(0)y(t) = DI > 0, lo que implica que xes creciente en x= 0. En conclusión, se tiene que las soluciones del sistema (6.1) son no negativas para cualquier t≥t0.
120 CAPÍTULO 6. MODELO SIN CRECIMIENTO EN PARED. Teorema 6.3.1 Supongamos que δ:= D−c > 0,(6.5) entonces para todo t > t0, las soluciones del sistema (6.1) están acotadas para cualquier condición inicial1. Demostración.- Definimos el funcional u(t) := u(x(t), y(t)) = x(t) + y(t). Derivando ua lo largo de soluciones del sistema (6.1) du(t) dt =dx(t) dt +dy(t) dt =DI −Dx(t)−aU(x(t))y(t)−Dy(t) + cU(x(t−ρ(t)))y(t) ≤DI −Dx(t)−ay(t)−Dy(t) + cy(t) ≤DI −Dx(t)−Dy(t) + cy(t) =DI −Dx(t)−(D−c)y(t) =DI −Dx(t)−δy(t), donde 0< δ < D definido como en (6.5). Definimos Ω := {(x, y)∈R2 +:Dx +δy ≤DI}. Vamos a estudiar cotas superiores de las soluciones caso por caso. (a) Si una trayectoria comienza en t=t0en R2 +\Ω, entonces el función u(x, y)a lo largo de las trayectorias que comienzan en este punto deben ser decrecientes para todo t≥t0tal que (x, y)∈R2 +\Ω, pues en tal zona, Dx +δy > DI, con lo que du(t) dt =DI −Dx(t)−δy(t)<0, es decir, udecrece. De esta forma, u(t) := u(x(t), y(t)) ≤u(x(t0), y(t0)) = x(t0) + y(t0)≤ |φ1|+|φ2|. 1Notemos que tal condición inicial siempre está acotada pues φes continua en [−h, 0].
6.4. ANÁLISIS DE LA ESTABILIDAD. 121 Así, k(x(t), y(t))k:= x(t) + y(t) = u(t)≤ |φ1|+|φ2|,(6.6) con lo que xeyestán acotadas. (b) Si una trayectoria comienza o entra en la región Ωen un tiempo t=t1≥t0y está allí para siempre, entonces por la definición de Ωtenemos que para t≥t1,Dx(t) + δy(t)≤DI, de donde se deduce D δx(t) + y(t)≤DI δ. De esta forma, como 0< δ < D, entonces D/δ > 1, con lo que k(x(t), y(t))k:= x(t) + y(t)<D δx(t) + y(t)≤DI δ,(6.7) de donde se tiene que xeyestán acotadas. (c) Si una trayectoria comienza y entra o vuelve a entrar después de haber salido de la región Ωen un tiempo t=t2i−1≥t0y sale en el tiempo t=t2i,(i= 1,2,...), entonces (6.6) se verifica para todos los tiempos (t2i, t2i+1)y (6.7) para todos los tiempos (t2i−1, t2i). En resumen, para cualquier t > t0, las soluciones (x(t), y(t)) están acotadas por k(x(t), y(t))k:= x(t) + y(t)≤m´ax |φ1|+|φ2|,DI δ. 6.4. Análisis de la estabilidad. Consideramos el modelo del quimiostato con retardo dado por el sistema (6.1) con una función de retardo variable ρ(·)∈C1(R; [0, h]). Es fácil darse cuenta de que (I, 0) es un equilibrio factible y es asintóticamente estable cuando D≥c. De hecho, cuando D≥c, se sigue directamente de la segunda ecuación del sistema (6.1) que y′(t)<0para todo t. Por tanto, y= 0 es un equilibrio asintóticamente estable de la segunda ecuación de (6.1), y haciendo y= 0 en la primera ecuación del sistema (6.1) tenemos inmediatamente x=I.
122 CAPÍTULO 6. MODELO SIN CRECIMIENTO EN PARED. Cualquier otro equilibrio no trivial2(x∗, y∗)debe satisfacer DI −Dx∗−ax∗ m+x∗y∗= 0,(6.8) −Dy∗+cx∗ m+x∗y∗= 0.(6.9) De aquí, dividiendo (6.9) por y∗y despejando x∗, tenemos cx∗ m+x∗=D⇔cx∗=mD +Dx∗⇔cx∗−Dx∗=mD ⇔ ⇔x∗(c−D) = mD ⇔x∗=mD c−D. Por otro lado, de (6.8), obtenemos DI −Dx∗−ax∗ m+x∗y∗= 0, es decir, DI −DmD c−D−a mD c−D m+mD c−D y∗= 0, de donde DI −mD2 c−D=a ✚ mD ✘✘ c−D ✚ mc−✟ ✟ mD+✟ ✟ mD ✘✘ c−D y∗, o lo que es lo mismo, DI −mD2 c−D=aD cy∗, o de otra forma y∗=DI −mD2 c−D aD c =cDI −cmD2 c−D aD = c2DI−cD2I−cmD2 c−D aD =c2DI −cD2I−cmD2 (c−D)aD . En resumen, x∗=mD c−D, y∗=c2DI −cD2I−cmD2 (c−D)aD .(6.10) 2Ahora sólo estamos interesados en los equilibrios positivos, debido al significado biológico de las soluciones de nuestro problema.
6.4. ANÁLISIS DE LA ESTABILIDAD. 123 Nos disponemos a analizar las propiedades de estabilidad del equilibrio positivo (x∗, y∗). Para ello construiremos un funcional de Lyapunov adecuado y obtendremos condiciones suficientes sobre los parámetros asociados a la función de retardo y sobre el resto de parámetros del modelo, con objeto de asegurar la estabilidad de este punto estacionario. Teorema 6.4.1 Supongamos que ρ(·) : R−→ [0, h]tal que 1−ρ′(t)≤M < 1. Entonces el equilibrio positivo (x∗, y∗)definido en (6.10) para el sistema (6.1) es globalmente asintóticamente estable. Nota importante.- Antes de realizar la demostración de este resultado, queremos hacer notar que, gracias al uso de técinas más sofisticadas tales como los funcionales de Lyapunov, hemos logrado mejorar el resultado análogo que aparece en [1], libro que estamos siguiendo a la hora estudiar esta parte del trabajo. Esto se debe a que en dicho libro se exige una cota sobre la amplitud del retardo, mientras que nosotros hemos logrado demostrar el mismo resultado sin pedir a cambio condición alguna. Demostración.- Primero, introducimos la siguiente transformación ex(t) := x(t)−x∗ ey(t) := log y(t) y∗(6.11) e U(·) := U(·+x∗)−U(x∗). Entonces, el sistema (6.1) puede escribirse como dex(t) dt =−Dex(t)−ay∗eey(t)e U(ex(t)) −ay∗U(x∗)eey(t)−1, dey(t) dt =ce U(ex(t−ρ(t))).(6.12) En efecto, como de (6.8) se tiene DI −Dx∗=ax∗ m+x∗y∗=aU(x∗)y∗, entonces dex(t) dt =DI −Dex(t)−Dx∗−ae U(ex(t)) + U(x∗)y∗eey(t) =−Dex(t)−ay∗eey(t)e U(ex(t)) −ay∗eey(t)U(x∗) + aU(x∗)y∗ =−Dex(t)−ay∗eey(t)e U(ex(t)) −ay∗U(x∗)eey(t)−1.
124 CAPÍTULO 6. MODELO SIN CRECIMIENTO EN PARED. Por otro lado, de (6.9) se deduce −Dy∗=−cx∗ m+x∗y∗=−cU(x∗)y∗, con lo que y∗eey(t)dey(t) dt =−Dy∗eey(t)+ce U(ex(t−ρ(t))) + U(x∗)y∗eey(t)(6.13) =−Dy∗eey(t)+ce U(ex(t−ρ(t)))y∗eey(t)+cU(x∗)y∗eey(t) =✭✭✭✭✭✭✭ ✭ −cU(x∗)y∗eey(t)+ce U(ex(t−ρ(t)))y∗eey(t)+✘✘✘✘✘✘ ✘ cU(x∗)y∗eey(t), y por tanto dey dt =−e U(ex(t−ρ(t))). Ahora construímos un funcional de Lyapunov adecuado paso a paso. Primero definimos V1(ex(t)) := Zex(t) 0e U(s)ds. Ahora bien, si x(t)> x∗, entonces ex(t)>0; de esta forma, como s∈(0,ex(t)), entonces s > 0. En consecuencia, s+x∗> x∗yUes una función creciente respecto de x, por lo que e U(s) = U(s+x∗)−U(x∗)>0. En conclusión, V1(ex(t)) = Zex(t) 0e U(s)ds > 0, siempre y cuando x(t)6=x∗. Análogamente, si x(t)< x∗, entonces ex(t)<0; por tanto, como s∈(ex(t),0), entonces s < 0. Entonces s+x∗< x∗yUes una función creciente respecto de x, con lo que e U(s) = U(s+x∗)− U(x∗)<0. En conclusión, V1(ex(t)) = Zex(t) 0e U(s)ds =−Zex(t)e U(s)ds > 0, siempre y cuando x(t)6=x∗. Resumiendo, V1(ex)>0para todo x6=x∗.
6.4. ANÁLISIS DE LA ESTABILIDAD. 125 Derivando V1respecto de ta lo largo de soluciones del sistema (6.12) tenemos dV1(t) dt =dex(t) dt e U(ex(t)) =−Dex(t)e U(ex(t)) −ay∗eey(t)e U2(ex(t)) −ay∗U(x∗)eey(t)−1e U(ex(t)). Definimos ahora V2(ey(t)) := Zey(t) 0 (es−1)ds. Ahora bien, si y(t)> y∗, entonces ey(t)>0; de esta forma, como s∈(0,ey(t)), entonces s > 0. En consecuencia, es−1>1−1 = 0, por lo que V2(ey(t)) = Zey(t) 0 (es−1) ds > 0, siempre y cuando y(t)6=y∗. Análogamente, si y(t)< y∗, entonces ey(t)<0; por tanto, como s∈(ey(t),0), entonces s < 0. Entonces es−1<1−1 = 0, con lo que V2(ey(t)) = Zey(t) 0 (es−1)ds =−Zey(t) (es−1)ds > 0, siempre y cuando y(t)6=y∗. Resumiendo, V2(ey)>0para todo y6=y∗. Derivando V2respecto de ta lo largo de soluciones del sistema (6.12) tenemos dV2(t) dt =dey(t) dt eey(t)−1 =ceey(t)−1e U(ex(t−ρ(t))). Consideramos ahora el funcional V(t) := V(ex(t),ey(t)) = V1(ex(t−ρ(t))) + λV2(ey(t)), con λ:= a c(1 −M)U(x∗)y∗. Está claro que V(ex(t),ey(t)) >0para todo (x, y)6= (x∗, y∗). Notemos que si (x, y) = (x∗, y∗), me encuentro en el equilibrio, por tanto, por definición de punto estacionario, no me moveré de él.
132 CAPÍTULO 7. MODELO CON CRECIMIENTO EN PARED. todos los tiempos (t2i, t2i+1)y (7.6) para todos los tiempos (t2i−1, t2i). En conclusión, para cualquier tiempo t > 0, tenemos que las soluciones de (x(t), y1(t), y2(t)) están acotadas por k(x(t), y1(t), y2(t))k:= x(t) + y1(t) + y2(t) <m´ax |φ1| b+1 + νh 1−M1|φ21|+|φ22|,DI bµ . 7.4. Análisis de la estabilidad. Notemos que cualquier solución de equilibrio (x, y1, y2)del sistema (7.1) debe satisfacer las ecuaciones algebraicas DI −Dx −aU(x)(y1+y2) + bνy1= 0, −(ν+D)y1+cU(x)y1−r1y1+r2y2= 0, −νy2+cU(x)y2+r1y1−r2y2= 0. Claramente (x, y1, y2) = (I, 0,0) es un equilibrio asintóticamente estable cuando ν≥c; de hecho, cuando ν≥c, d(y1(t) + y2(t)) dt =−(ν+D)y1(t)−νy2(t) + cU(x(t−τ2(t))(y1(t) + y2(t)) <−(ν−c)(y1(t) + y2(t)) <0. Para estudiar la estabilidad de las soluciones de equilibrio no triviales del sistema (7.1), primero llevamos a cabo el siguiente cambio de variables: η(t) := y1(t) + y2(t); β(t) := y1(t) η(t). De esta forma, y análogamente a lo estudiado en los primeros capítulos de esta memoria, el sistema (7.1) se reescribe como dx(t) dt =DI −Dx(t)−aU(x(t))η(t) + bνβ(t−τ1(t))η(t−τ1(t)), dη(t) dt =−νη(t)−Dβ(t)η(t) + cU(x(t−τ2(t)))η(t),(7.7) dβ(t) dt =−Dβ(t)(1 −β(t)) −r1β(t) + r2(1 −β(t)).
7.4. ANÁLISIS DE LA ESTABILIDAD. 133 La ecuación para β(t)está desacoplada de las ecuaciones para x(t)yη(t). Resolviendo la ecuación algebraica −Dβ(1 −β)−r1β+r2(1 −β) = 0 tenemos, como ya sabemos, β=D+r1+r2±p(D+r1+r2)2−4Dr2 2D. Notemos que β(t)∈(0,1) para cualesquiera soluciones positivas y1(t)>0,y2(t)>0del sistema (7.1). Además, como β′β=0 =r2>0yβ′β=1 =−r1<0, entonces (0,1) es positivamente invariante para β(t)y es la región biológicamente interesante1. Por tanto, el equilibrio de la ecuación desacoplada del sistema (7.6) está dado por β∗:= D+r1+r2−p(D+r1+r2)2−4Dr2 2D.(7.8) Lema 7.4.1 El punto crítico β∗tal como está definido en (7.8) es un equilibrio asintóticamente estable de la tercera ecuación del sistema (7.6). Demostración.- Definimos ∆(t) := β(t)−β∗, entonces ∆(t)satisface la siguiente EDO: d∆(t) dt =−D(∆(t) + β∗)(1 −∆(t)−β∗)−r1(∆(t) + β∗) + r2(1 −∆(t)−β∗) =−D∆(t)(1 −∆(t)−β∗)−Dβ∗(1 −∆(t)−β∗)−r1∆(t)−r1β∗ +r2−r2∆(t)−r2β∗ =−D∆(t) + D∆2(t) + D∆(t)β∗−Dβ∗+Dβ∗∆(t) +Dβ∗2−r1∆(t)−r1β∗+r2−r2∆(t)−r2β∗. Ahora bien, por ser β∗equilibrio de la tercera ecuación del sistema (7.6), se tiene −Dβ∗(1 −β∗)−r1β∗+r2(1 −β∗) = 0, es decir, −Dβ∗+Dβ∗2−r1β∗+r2−r2β∗= 0. Continuando d∆(t) dt =−D∆(t) + D∆2(t) + D∆(t)β∗+Dβ∗∆(t)−r1∆(t)−r2∆(t) =D∆(t)∆(t) + 2β∗−1−r1+r2 D, 1En este párrafo no damos más detalles por analogía a estudios anteriores
134 CAPÍTULO 7. MODELO CON CRECIMIENTO EN PARED. es decir, d∆(t) dt =D∆(t) ∆(t)−p(D+r1+r2)2−4Dr2 D!.(7.9) Resolviendo (7.9) con condición inicial ∆(t0) = ∆0, tenemos ∆(t) = ∆0p(D+r1+r2)2−4Dr2 ∆0−(∆0−p(D+r1+r2)2−4Dr2 e√(D+r1+r2)2−4Dr2(t−t0)−→ 0,(7.10) cuando t−→ 0, lo que implica que β∗es asintóticamente estable. Para finalizar, veamos que se tiene (7.10). En efecto, partimos de (7.9), que podemos reescribir como d∆(t) dt =D∆2(t)−p(D+r1+r2)2−4Dr2∆(t). Para un tiempo suficientemente grande, reemplazando β(t)por β∗en las dos primeras ecuaciones del sistema (7.6), tenemos dx(t) dt =DI −Dx(t)−aU(x(t))η(t) + bνβ∗η(t−τ1(t)), dη(t) dt =−νη(t)−Dβ∗η(t) + cU(x(t−τ2(t)))η(t).(7.11) Cualquier solución de equilibrio no trivial del sistema (7.11) debe satisfacer las ecuaciones algebraicas DI −Dx∗−aU(x∗)η∗+bνβ∗η∗= 0,(7.12) −ν−Dβ∗+cU(x∗==0.(7.13) Despejando x∗en (7.13), tenemos −ν−Dβ∗+cx∗ m+x∗= 0, es decir, cx∗ m+x∗=ν+Dβ∗,
7.4. ANÁLISIS DE LA ESTABILIDAD. 135 o lo que es lo mismo, cx∗=mν +mDβ∗+x∗ν+Dβ∗x∗. Reagrupando los términos con x∗llegamos a cx∗−νx∗−Dβ∗x∗=mν +mDβ∗; es decir, (c−ν−Dβ∗)x∗=m(ν+Dβ∗), con lo que x∗=m(ν+Dβ∗) c−ν−Dβ∗. Ahora, como sabemos que −ν−Dβ∗+cU(x∗) = 0, entonces cU(x∗) = ν+Dβ∗, con lo que U(x∗) = ν+Dβ∗ c y sustituyendo esta expresión en (7.12), tenemos DI −Dx∗−aν+Dβ∗ cη∗+bνβ∗η∗= 0, es decir, DI −Dx∗−ha c(ν+Dβ∗)−bνβ∗iη∗= 0, con lo que η∗=DI −Dx∗ a c(ν+Dβ∗)−bνβ∗=DI −Dx∗ a c(ν+Dβ∗)−bcνβ∗ c =cDI −cDx∗ a(ν+Dβ∗)−bcνβ∗=cD(I−x∗) a(ν+Dβ∗)−bcνβ∗.
136 CAPÍTULO 7. MODELO CON CRECIMIENTO EN PARED. Ahora bien, para que sea equilibrio positivo, debe ocurrir que c−ν−Dβ∗>0, i.e., c > ν +Dβ∗ y también que cD(I−x∗)(a(ν+Dβ∗)−bcνβ∗)>0. Ahora bien, I−x∗=I−m(ν+Dβ∗) c−ν−Dβ∗=I(c−ν−Dβ∗)−m(ν+Dβ∗) c−ν−Dβ∗, por tanto, como c−ν−Dβ∗>0(es la primera condición exigida), entonces la segunda condición es equivalente a (I(c−ν−Dβ∗)−m(ν+Dβ∗))(a(ν+Dβ∗)−bcνβ∗)>0. A continuación establecemos un resultado de estabilidad para el equilibrio (x∗, η∗). Teorema 7.4.1 Supongamos que τ1(·), τ2(·) : R+−→ [0, h],τ′ i(t)≤Mi<1,i= 1,2; entonces el equilibrio (x∗, η∗)del sistema (7.11) es asintóticamente estable, siempre y cuando 1 ν+Dβ∗e(ν+Dβ∗)h−1+1 1−M1 h < 2a bcνβ∗. Demostración.- En primer lugar, consideramos la transformación ex(t) := x(t)−x∗, eη(t) := log η(t) η∗, e U(·) := U(·+x∗)−U(x∗). De esta manera, el sistema (7.11) queda reescrito como dex(t) dt =−Dex(t)−ae U(ex(t))η∗eeη(t)−aU(x∗)η∗(eeη(t)−1) + bνβ∗η∗(eeη(t−τ1(t)) −1), deη(t) dt =ce U(ex(t−τ2(t))). Veamos que, efectivamente, es cierto. Antes de nada, de (7.12) tenemos DI −Dx∗=aU(x∗)η∗−bνβ∗η∗.
7.4. ANÁLISIS DE LA ESTABILIDAD. 137 Entonces, dex(t) dt =DI −Dex(t)−Dx∗−a(e U(ex(t)) + U(x∗))η∗eeη(t)+bνβ∗η∗eeη(t−τ1(t)) =DI −Dex(t)−Dx∗−ae U(ex(t))η∗eeη(t)−aU(x∗)η∗eeη(t)+bνβ∗η∗eeη(t−τ1(t)) =aU(x∗)η∗−bνβ∗η∗−Dex(t)−ae U(ex(t))η∗eeη(t)−aU(x∗)η∗eeη(t)+bνβ∗η∗eeη(t−τ1(t)) =−Dex(t)−ae U(ex(t))η∗eeη(t)−aU(x∗)η∗eeη(t)−1+bνβ∗η∗eeη(t−τ1(t)) −1. Por otro lado, de (7.13) se tiene −ν−Dβ∗+cU(x∗) = 0, entonces η∗eeη(t)deη(t) dt =−νη∗eeη(t)−Dβ∗η∗eeη(t)+c(e U(ex(t−τ2(t))) + U(x∗))η∗eeη(t) =−νη∗eeη(t)−Dβ∗η∗eeη(t)+ce U(ex(t−τ2(t)))η∗eeη(t)+cU(x∗)η∗eeη(t) = (−ν−Dβ∗+cU(x∗))η∗eeη(t)+ce U(ex(t−τ2(t)))η∗eeη(t) =ce U(ex(t−τ2(t)))η∗eeη(t). Así, deη(t) dt =ce U(ex(t−τ2(t))). Vamos a definir de nuevo un funcional de Lyapunov adecuado paso a paso. Definimos V1(ex(t)) := Zex(t) 0e U(s)ds. Ahora bien, si x(t)> x∗, entonces ex(t)>0; de esta forma, como s∈(0,ex(t)), entonces s > 0. En consecuencia, s+x∗> x∗yUes una función creciente respecto de x, por lo que e U(s) = U(s+x∗)−U(x∗)>0. En conclusión, V1(ex(t)) = Zex(t) 0e U(s)ds > 0, siempre y cuando x(t)6=x∗. Análogamente, si x(t)< x∗, entonces ex(t)<0; por tanto, como s∈(ex(t),0), entonces s < 0. Entonces s+x∗< x∗yUes una función creciente respecto de x, con lo que e U(s) = U(s+x∗)− U(x∗)<0. En conclusión,
138 CAPÍTULO 7. MODELO CON CRECIMIENTO EN PARED. V1(ex(t)) = Zex(t) 0e U(s)ds =−Zex(t)e U(s)ds > 0, siempre y cuando x(t)6=x∗. Resumiendo, V1(ex)>0para todo x6=x∗. Derivando V1respecto de ta lo largo de soluciones del sistema (7.11) tenemos dV1(t) dt =dex(t) dt e U(ex(t)) =−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−aU(x∗)η∗eeη(t)−1e U(ex(t)) +bνβ∗η∗eeη(t−τ1(t)) −1e U(ex(t)) =−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−aU(x∗)η∗eeη(t)−1e U(ex(t)) +bνβ∗η∗eeη(t−τ1(t)) −1e U(ex(t)) ±bνη∗β∗eeη(t)−1e U(ex(t)) =−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) −bνβ∗η∗eeη(t)−eeη(t−τ1(t))e U(ex(t)) −✘✘✘✘✘✘✘ ✘ bνβ∗η∗e U(ex(t)) + ✘✘✘✘✘✘✘ ✘ bνβ∗η∗e U(ex(t)) =−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) −bνβ∗η∗eeη(t)−eeη(t−τ1(t))e U(ex(t)) <−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) +bνβ∗η∗eeη(t)−eeη(t−τ1(t))e U(ex(t)). Ahora bien, eeη(t)−eeη(t−τ1(t))=Zt t−τ1(t) d ds eeη(s)ds=Zt t−τ1(t) eeη(s)deη(s) ds ds =cZt t−τ1(t) eeη(s)e U(ex(s−τ2(s)))ds≤cZt t−τ1(t)eeη(s)e U(ex(s−τ2(s)))ds =cZt t−τ1(t) eeη(s)e U(ex(s−τ2(s)))ds y utilizando la desigualdad de Young
7.4. ANÁLISIS DE LA ESTABILIDAD. 139 dV1(t) dt <−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) +bνβ∗η∗eeη(t)−eeη(t−τ1(t))e U(ex(t)) ≤ −Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) +bcνβ∗η∗Zt t−τ1(t) eeη(s)e U(ex(s−τ2(s)))dse U(ex(t)) =−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) +bcνβ∗η∗Zt t−τ1(t) eeη(s)e U(ex(s−τ2(s)))e U(ex(t))ds ≤ −Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) +bcνβ∗η∗ 2Zt t−τ1(t) eeη(s)ds ·e U2(ex(t)) + Zt t−τ1(t) eeη(s)e U(ex(s−τ2(s))) 2ds =−Dex(t)e U(ex(t)) −ae U2(ex(t))η∗eeη(t)−(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t)) +bcνβ∗η∗ 2Zt t−τ1(t) eeη(s)ds ·e U2(ex(t)) + Zt t−τ1(t) eeη(s)e U2(ex(s−τ2(s)))ds En segundo lugar definimos V2(eη(t)) := Zeη(t) 0 (es−1)ds. Ahora bien, si η(t)> η∗, entonces eη(t)>0; de esta forma, como s∈(0,eη(t)), entonces s > 0. En consecuencia, es−1>1−1 = 0, por lo que V2(eη(t)) = Zeη(t) 0 (es−1) ds > 0, siempre y cuando η(t)6=η∗. Análogamente, si η(t)< η∗, entonces eη(t)<0; por tanto, como s∈(eη(t),0), entonces s < 0. Entonces es−1<1−1 = 0, con lo que V2(eη(t)) = Zeη(t) 0 (es−1)ds =−Zeη(t) (es−1)ds > 0, siempre y cuando η(t)6=η∗. Resumiendo, V2(eη)>0para todo η6=η∗.
140 CAPÍTULO 7. MODELO CON CRECIMIENTO EN PARED. Derivando V2respecto de ta lo largo de soluciones del sistema (7.11) tenemos dV2(t) dt =deη(t)) dt eeη(t)−1e U(ex(t−τ2(t))). Finalmente, definimos el funcional V3(ex(t),eη(t)) := 1−M2 2(1 −M1)bcνη∗β∗Zt t−τ1(t)Zt v eeη(s)e U2(ex(s−τ2(s)))dsdv. Claramente V3(ex, eη)>0para todo η6=η∗,x6=x∗, ya que todas las constantes que aparecen, así como el integrando, son estrictamente positivos, teniendo también en cuenta los extremos de integración. Derivando V3respecto de ta lo largo de soluciones del sistema (7.11) tenemos dV3(t) dt =1−M2 2(1 −M1)bcνη∗β∗d dt Zt t−τ1(t)Zt ve U2(ex(s−τ2(s)))dsdv =1−M2 2(1 −M1)bcνη∗β∗Zt t−τ1(t)d dt Zt v eeη(s)e U2(ex(s−τ2(s)))dsdv −Zt t−τ1(t) eeη(s)e U2(ex(s−τ2(s)))·(1 −τ′ 1(t)) =1−M2 2(1 −M1)bcνη∗β∗neeη(s)e U2(ex(s−τ2(s)))τ1(t) −Zt t−τ1(t) eeη(s)e U2(ex(s−τ2(s)))ds·(1 −τ′ 1(t)). Ahora consideramos el funcional V(t) = V(ex(t),eη(t)) = V1(ex(t−τ2(t))) + κV2(eη(t)) + V3(eη(t)), con κ=(aU(x∗)−bνβ∗)η∗ c(1 −M2). Derivando Vrespecto de ta lo largo de soluciones del sistema (7.11) y usando que τ′ i(t)≤ M1<1,i= 1,2, tenemos
7.4. ANÁLISIS DE LA ESTABILIDAD. 141 dV (t) dt =dV1(ex(t−τ2(t))) dt +κdV2(eη(t)) dt +dV3(eη(t)) dt =dex(t−τ2(t)) dt e U(ex(t−τ2(t))) + κdV2(eη(t)) dt +dV3(eη(t)) dt =dex(s) ds s=t−τ2(t) d(t−τ2(t)) dt e U(ex(t−τ2(t))) + κdV2(eη(t)) dt +dV3(eη(t)) dt =dex(s) ds s=t−τ2(t)·(1 −τ′ 2(t))e U(ex(t−τ2(t))) + κdV2(eη(t)) dt +dV3(eη(t)) dt =dV1(ex(s)) ds s=t−τ2(t)·(1 −τ′ 2(t)) + κdV2(eη(t)) dt +dV3(eη(t)) dt ≤n−Dex(t−τ2(t))e U(ex(t−τ2(t))) −ae U2(ex(t−τ2(t)))η∗eeη(t) −(aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t−τ2(t))) +bcνβ∗η∗ 2Zt t−τ1(t) eeη(s)ds ·e U2(ex(t−τ2(t)) + Zt t−τ1(t) eeη(2) e U2(ex(s−τ2(s)))ds· ·(1 −τ′ 2(t)) + (1 −M2)(aU(x∗)−bνβ∗)η∗ cceeη(t)−1e U(ex(t−τ2(t))) +1−M2 2(1 −M1)bcνη∗β∗neeη(t)e U2(ex(t−τ2(t)))τ1(t) −Zt t−τ1(t) eeη(s)e U2(ex(s−τ2(s)))ds·(1 −τ′ 1(t)) ≤n−Dex(t−τ2(t))e U(ex(t−τ2(t))) −ae U2(ex(t−τ2(t))η∗eeη(t) −✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ (aU(x∗)−bνβ∗)η∗eeη(t)−1e U(ex(t−τ2(t))) +bcνβ∗η∗ 2Zt−τ1(t) eeη(s)ds ·e U2(ex(t−τ2(t))) +✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ ✭ Zt−τ1(t) eeη(s)e U2(ex(s−τ2(s)))ds· ·(1 −M2) + ✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ ✭ (1 −M2)(aU(x∗)−bνβ∗)η∗ cceeη(t)−1e U(ex(t−τ2(s))) +1−M2 2(1 −M1)bcνη∗β∗neeη(t)e U2(ex(t−τ2(t)))h −✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ ✭ (1 −M1)Zt t−τ1(t) eeη(s)e U2(ex(s−τ2(s)))ds) =−Dex(t−τ2(t))e U(ex(t−τ2(t))) −ae U2(ex(t−τ2(t)))η∗eeη(t) +bcνη∗β∗ 2Zt t−τ1(t) eeη(s)ds ·e U2(ex(t−τ2(t))) +1−M2 2(1 −M1)bcνη∗β∗eeη(t)e U2(ex(t−τ2(t)))h. (7.14)