Continuous Dynamical Systems Florence Hubert - - PowerPoint PPT Presentation

Continuous Dynamical Systems

Florence Hubert

florence.hubert@univ-amu.fr

Charlotte Perrin

charlotte.perrin@univ-amu.fr

Master CMB-BIO,

2019-2020

Practical information

schedule

◮ lectures + exercises 12h ◮ computer sessions (language Python) 2 × 3h

evaluation

◮ group project + reports on the computer sessions ◮ final exam

Mathematics and applications of mathematics

http://mathsmonde.math.cnrs.fr https://images.math.cnrs.fr/

Basics in the modeling methodology

Example of simple biological system: population of individuals unknown and variable

◮ the size of the population N is seen as a function of time: N = N(t)

discrete description of the evolution of the system

◮ observation at discrete times t0 = 0, t1 = t0 + ∆t, t2 = t0 + 2∆t, etc.

∆t = constant time interval separating two observations

◮ variation of the population between two observations is given by

N(ti+1) − N(ti) = nb of births − nb of deaths ← conservation principle

◮ we prescribe how the nb of deaths and births (during ∆t) is related to N

nb of births = λ ∆t N(ti), nb of deaths = µ ∆t N(ti) ← comportemental law (Malthus)

passage to the continuous description ∆t → 0 d dt N(t) = aN(t) where a = λ − µ

Continuous dyn. systems Differential Equations (ODEs)

More generally ... (1st order) ODE problem: given a function f , find t → y(t) satisfying d dt y(t) = f

• t, y(t)
• ◮ when f only depends on y, the ODE is said to be autonomous

◮ when f depends linearly on y, the ODE is said to be linear

Cauchy problem: given t0, y0 (and f ), find y satisfying        d dt y(t) = f

• t, y(t)
• y(t0) = y0

◮ y0 is called the initial condition / initial datum

“Well-posed” problem

• J. Hadamard (1865-1963)

→ consistency of the models

1

existence of a solution

2

uniqueness of the solution

3

continuous dependency with respect to the (initial) data (/ stability) well-posedness of the Cauchy problem → Cauchy-Lipschitz Theorem

Notion of solution

Cauchy problem: given t0 ∈ I ⊂ R, y0 ∈ U ⊂ Rn, find y solution to        d dt y(t) = f

• t, y(t)
• y(t0) = y0

Definition (Notion of solution)

A solution to the Cauchy problem is a couple (J, y) where J ⊂ I and y : J → U is a C1 function satisfying the ODE for all t ∈ J.

Definition (Global solution)

A solution to the Cauchy problem is called global if J = I.

Example: solution to the Malthus model

       d dt N(t) = aN(t) N(0) = N0 Linear, autonomous (1st order) ODE

Example: solution to the Malthus model

       d dt N(t) = a N(t) N(0) = N0 Linear, autonomous (1st order) ODE global solution (R, N) given by

N(t) = N0eat ∀t ∈ R

Representation of the solution

Cauchy problem: given t0 ∈ I ⊂ R, y0 ∈ U ⊂ Rn, find y solution to        d dt y(t) = f

• t, y(t)
• y(t0) = y0

Definition (Notion of solution curves)

Let (J, Y ) a solution of the problem, we call solution curve of the system the set {(t, Y (t)), t ∈ J} ⊂ R × Rn.

Definition (Notion of trajectories or orbits)

Let (J, Y ) a solution of the problem, we call trajectory or orbit of the system the set {Y (t), t ∈ J} ⊂ Rn.

Example: prey predator model

• x′(t) = f1
• x(t), y(t)
• y ′(t) = f2
• x(t), y(t)
• x(0) = x0

y(0) = y0 number of preys x and predators y for different initial data left: the solution curves (x in plain lines, y in dotted lines), right: the trajectories

I II III IV

Some mathematical issues addressed in this course

local (and global) well-posedness equilibrium points y ∗ which are such that f (t, y ∗) = 0 for all t stability of equilibrium points long-time behavior

• ther interesting issues:

numerical discretization, dependency with respect to parameters and bifurcation phenomena, etc.

Outline

ODEs and biology: most famous models linear systems of ODEs global well-posedness, stability of the origin, long-time behavior nonlinear ODEs conditions for local well-posedness, possible blow-up phenomena, equilibrium points and their stability

Some ODEs arising in biology

Population dynamics Epidemiology Biochemistry Electrical networks in biology

Population dynamics: Malthusian growth model

       d dt N(t) = a N(t) N(0) = N0 Linear, autonomous (1st order) ODE global solution (R, N) given by

N(t) = N0eat ∀t ∈ R

exponential growth or decrease depending on the sign of a ref: Malthus (1798)

Malthusian catastrophe

1944: introduction of 29 reindeer on St Matthew island (Alaska) 1963: more than 6000 individuals observed (increase of more than 30% per year) 1966: almost the entire population is dead of starvation. ❀ Lichens had been completely eliminated ! ! the model does not take into account the possible limitation of the ressources ! A Malthusian catastrophe refers to a demographic collapse that follows an exponential growth of the population.

Models for one population with self-limiting growth

→ The Logistic model

Malthus: d dt N(t) = a N(t)

• d

dt N(t) = a

• 1 − N(t)

b

• N(t)

b is the carrying capacity which depends on the environment non-linear 1st order (autonomous if a and b are constant) ODE equilibrium points: N∗

1 = 0, N∗ 2 = b

there exists a global solution and an analytical expression for N can be derived N(t) → b as t → +∞ sigmoid shape of the sol. curves ref: Verhulst (1838)

1 2 3 4 5 5 10 15 20 25 30 Logistique, a=1, b=10 Logistique, a=1, b=20 Logistique, a=2, b=10 Logistique, a=2, b=20

Models for one population with self-limiting growth

→ Gombertz and others..

Gombertz (1825): introduced for insurance companies, today in cancer modeling N′(t) = a ln

• b

N(t)

• N(t)

◮ analogue prop as the logistic model

2 4 6 8 10 5 10 15 20 25 30 Logistique, a=1, b=10 Logistique, a=1, b=20 Gompertz, a=1, b=10 Gompertz, a=1, b=20

more generally: N′(t) = h(N(t))N(t)

◮ h(N) = Nγ with γ < 0: power law model; ◮ h(N) = aN− 1 4 − b: West’s model; ◮ h(N) = aN− 1 3 − b: Von Bertalanffy’s model.

Models for several populations: the prey-predator model

Biological assumptions in the absence of predators, the preys have an unlimited growth; the survival of the predators depends on the preys; the growth of the predator population is proportional to the size of the prey population. x: number of preys, y: number of predators, parameters a, b, c, d ≥ 0

• x′(t) = a x(t) − b x(t)y(t)

y ′(t) = d x(t)y(t) − c y(t)

preys (water lilies) - predators (ducks)

• x′(t) = a x(t) − b x(t)y(t)

y ′(t) = d x(t)y(t) − c y(t)

autonomous (the parameters are cst) system of two coupled non-linear ODEs no explicit solution

I II III IV Left: solution curves, Right: corresponding trajectories in the xy-plane.

Extensions to more general predation rates

• x′(t) = ax(t) − bx(t)y(t)

y ′(t) = φ(x(t))y(t) − cy(t)

Density of prey population

Number of prey consumed I II III

Holling Type I: φ(x) = d

• x1{x<¯

x} + ¯

x1{x≥¯

x}

• ◮ the time needed for the predators to find and consume preys is negligible

Holling Type II: φ(x) =

dx 1+hdx

◮ h > 0: handling time of consumption ◮ limited capacity of the predator to find and eat the preys

Holling Type III: φ(x) =

dx2 1+hdx2

◮ for large x, saturation effect similar to type II ◮ for small x, the predation rate remains low → preys can evade

Two populations in competition

Two populations in the same environment exploiting the same limited resources Biological assumptions each population has a negative influence on the growth rate of the other each population satisfies a logistic growth in the absence of the other        x′

1(t) = a1 x1(t)

• 1 − x1(t)

b1 −αx2(t) b1

• x′

2(t) = a2 x2(t)

• 1 − x2(t)

b2 −β x1(t) b2

• a1,2 ≥ 0 are the intrinsic growth rates

b1,2 ≥ 0 are the carrying capacities α, β > 0 (resp. < 0) intensity of the competition (resp. symbiosis)

A model in epidemiology: the SIR model

Time evolution of a population split into three parts

population S of susceptible individuals; population I of infected individuals; population R of recovered (or immune) individuals.

Biological assumptions

the total population is constant; the average number of new infected individuals par unit time is proportional to the product of the infected and susceptible population sizes.      S′(t) = −βS(t)I(t) + γR(t) I ′(t) = βS(t)I(t) − δI(t) R′(t) = δI(t) − γR(t) β: infection rate δ: recovery rate γ: rate of transition between R and S

Basic reproduction number R0

If R0 = β δ N > 1, the infection will be able to spread in the population. prediction of the spread of a disease possible extension to take into account the effect of a vaccination campaign...

Pharmacokinetics (PK)/Pharmacodynamics (PD) of a drug

1

PK : How the organism affects the drug

2

PD : How the drug affects the organism

Pharmacokinetics of a drug

How the organism affects the drug

One compartment model - infusion- one administration

c1 = Q1

V1 : concentration

k10 = Cl1

V1 : elimination rate

   dc1 dt (t) = −Cl1 V1 c1(t) + D(t) V1 c(tinj) = 0 t → D(t): injection protocol Cl: Clearance

Pharmacokinetics of a drug

One compartment model - Oral administration - one tablet

Example of targeted therapy (kinase inhibitors): Imatinib

◮ Chronic myelogenous leukemia or Gastro Intestinal Stromal Tumors

dqa dt = −kaqa, qa(tabs) = D dc1 dt = − Cl V1 c + ka V1 qa, c1(tabs) = 0 qa: quantity of absorbed drug

Pharmacokinetics of a drug

Two compartments model - injection Example of anti-angiogenic drug : Bevacizumab

ref: Bruno et al 1996

◮ Lung cancer, kidney cancer, glioblastoma,...

dc1 dt = − (k10 + k12) c1 + k21 V2 V1 c2(t) + u(t) V1 dc2 dt = k12 V1 V2 c1 − (k20 + k21) c2

Chemical kinetics

Elementary reactions α1A1 + α2A2 + · · · + αpAp

k

− → β1B1 + β2B2 + · · · + βrBr unknowns: the reactant concentrations [Ai], the product concentrations [Bi] parameters: αi, βi the stoichiometric coefficients, k the kinetic constant ODE system d dt [Ai] = −αik

p

• j=1

[Aj]αj and d dt [Bi] = βik

p

• j=1

[Aj]αj.

Chemical kinetics

The simplest reaction A

k

− → B ❀ Linear 1st order equation d dt [A] = −k[A] = − d dt [B] Order 1 law = ⇒ [A](t) = [A](0)e−kt → exponential decrease The case of two reactants and one product A + B

k

− → C ❀ Non-linear 1st order differential system d dt [C] = k[A][B] = − d dt [A] = − d dt [B] Note that d dt [A] = k[A](Cst − [A])

Complex chemical reaction

Example 1 A

k1

− → B, B

k2

− → C Each reaction can be represented by an order 1 law      [A]′(t) = −k1[A](t) [B]′(t) = k1[A](t) − k2[B](t) [C]′(t) = k2[B](t) ❀ Linear ODE system

Complex chemical reaction

Example 2: H2 + I2 − → 2HI The dynamics can be described by the elementary reactions I2 ⇋ I + I I + H2 ⇋ IH2 I + IH2 − → HI + HI We get then the differential system d dt [I2] = −k1[I2] + k−1[I]2 d dt [I] = 2k1[I2] − k−1[I]2 − k2[I][H2] + k−2[IH2] − k3[I][IH2] d dt [H2] = −k2[I][H2] + k−2[IH2] d dt [IH2] = k2[I][H2] − k−2[IH2] − k3[I][H2] d dt [HI] = 2k3[I][IH2]

Complex chemical reaction

Example 3: enzyme kinetics E + S

kon

⇋

koff ES kcat

→ E + P

a substrate S binds reversibly to an enzyme E to form an enzyme-substrate complex ES ES then reacts irreversibly to generate a product P and to regenerate the free enzyme E d dt [S] = −kon[E][S] + koff [ES] d dt [E] = −kon[E][S] + (koff + kcat)[ES] d dt [P] = kcat[ES] d dt [ES] = kon[E][S] − (koff + kcat)[ES] If [ES] reaches rapidly a steady state, we obtain the Michaelis-Menten equation d dt [S] = v = kcat[ET]

koff +kcat kon

+ [S] [S] = VM K + [S][S].

Gene circuit models - ref: Steinacher 2016

Understand how organisms can be at the same time robust to small perturbations (noise) and also sufficiently reactive to respond to evolutionary pressures coming from their environment. Biological processes of the gene expression the transcription: this is the production of the RNA copy of the DNA. This process creates messenger RNA (mRNA); the translation: a mRNA molecule (the mRNA concentration is denoted M), is translated into many protein molecules (the protein concentration is denoted P); the mRNA and protein degradation. Two models for the activation of the transcription process

Auto-activation

                       M′(t) = a  k1 + k2 P(t) kD N   1 + P(t) kD N − k4M(t) P′(t) = k3M(t) − k5P(t)

Simple activation

                       M′(t) = a  k1 + k2 T kD N   1 + T kD N − k4M(t) P′(t) = k3M(t) − k5P(t)

Hill function H(x) = 1 + 2 x kD N 1 + x kD N

a maximum transcription rate, k1 basal transcription rate, k2 the rate of feedback-mediated transcription, k3 translation rate, k4 and k5 mRNA and protein degradation rates, N Hill coefficient, kD link to half life of P.

❀ Asymptotic behaviour: mono or bistability in the auto-activated model, and

monostability in the simple activation model

Microtubules and RLC circuit

MTs are cytosketal polymers involved in many biological processes (cell division, cell migration, intracellular transport,...) Biological assumptions MTs seen as electric transmission lines with capacitive, resistive and inductive characteristics The capacitive and resistive components are due to the depleted layer which separates the interior of the MT from the bulk solution. The electric current propagating in a MT generates a magnetic flux inducing a voltage which is modeled by the inductive component

RLC circuit associated in series

Rt effective resistance (measuring the dissipation energy in the circuit), L the inductance of the inductor, C the capacitance (capacity of storing electric charges) of the capacitor. Uc UR, Ul the voltage across the capacitor the resistance and the inductor respectively.. I the current E electric generator

RLC circuit associated in series

Kirchhoff’s voltage law:

E = Uc + UL + UR;

Law of the resistance:

UR = RtI;

Law of the inductor:

UL = LdI dt ;

Law of the capacitor:

I = C dUc dt .

❀ Second order linear differential equation LC d2Uc dt2 + RtC dUc dt + Uc = E

RLC circuit associated in series

Y ′ = AY +b with A = 1 − 1

LC

− Rt

L

• and

b =

• E
• .

Case E = 0 (homogeneous system) and C = 1

The van der Pol oscillator

Observations in electrical circuits with vacuum tubes (triodes) Stable oscillations composed by rapid variations of the electrical potential, followed by a slower relaxation phase. ❀ relaxation-oscillation phenomenon An electric model for the heart Can be used to model the successive depolarization and repolarization phases (relaxation phase - slow process) of the heart cells membranes.

The van der Pol oscillator

Observations in electrical circuits with vacuum tubes (triodes) Stable oscillations composed by rapid variations of the electrical potential, followed by a slower relaxation phase. ❀ relaxation-oscillation phenomenon van der Pol system x′′(t) − 1 RC (1 − x2(t))x′(t) + 1 LC x(t) = 0

• r

     x′(t) = x(t) − 1 3x3(t) − y(t) y ′(t) = 1 T x(t) with T = L R2C . Oscillatory behaviour

In blue: T = 1 small, in red T = 10 large

FitzHugh-Nagumo system

Modeling the electric potential generated by a neuron      x′(t) = x(t) − 1 3x3(t) − y(t) + I(t) y ′(t) = 1 T (x(t) + a − by(t)) x(t) electrical potential of the axon membrane y(t) recovery variable I(t) stimulus current generated by the other neurons

FitzHugh-Nagumo system

Modeling the electric potential generated by a neuron      x′(t) = x(t) − 1 3x3(t) − y(t) + I(t) y ′(t) = 1 T (x(t) + a − by(t)) x(t) electrical potential of the axon membrane y(t) recovery variable I(t) stimulus current generated by the other neurons

In blue: T = 1 small, in red T = 10 large a = b = 1, I = 1.5