Formal Computational Modelling of Bone Physiology and Disease - - PowerPoint PPT Presentation

Formal Computational Modelling of Bone Physiology and Disease Processes

Candidate: Nicola Paoletti Supervisor: Emanuela Merelli Examiners: Luca Bortolussi, Guido Sanguinetti

PhD course in Information Science and Complex Systems - XVI cycle - University

• f Camerino

Camerino, 24 March 2014

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Outline

METHODS

Sect 1.1 Sect 1.2 Sect 2

D − CGF (extension

• f CGF PA)

Shape Calculus (Spatial PA) Controlled Switched Systems Stochastic Hybrid Automata Hybrid Ap- proximation Stochastic ABM CTMC ODE Simulation Probabilistic Model Checking Stabili- zation LANGUAGES MODELS ANALYSIS TECHNIQUES BIOLOGICAL PROPERTIES Simulation Sensitivity Analysis Model Checking Parameter Synthesis Optimal Control Bone micro structure Osteopo- rosis Osteomy- elitis STI (Therapy Scheduling)

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

The bone remodelling process

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

The bone remodelling process

Bone remodelling. . . why?

• Multiscale process: multiscale effects, from the molecular

signalling level to the tissue level

• A paradigm for other physiological systems (like

epithelium renewal and haematopoiesis)

• Clinical and social impact: for diseases like osteoporosis,

Paget’s disease, osteomyelitis.

• Collaboration with IOR (Istituto Ortopedico Rizzoli)

Shape Calculus (Spatial PA) Stochastic ABM Simulation Bone micro structure Osteopo- rosis

1.1- Shape Calculus & agent-based simulation

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Shape Calculus

• A bio-inspired spatial process algebra for describing 3D

processes, entities characterized by a behaviour B (` a la Timed CCS) and a shape S.

• We extend the calculus with additional terms: Iteration,

Thanatos and Duration

S ∈ S - 3D shape σ = V , m, p, v Basic shape SXS Complex shape B ∈ B - Behavior nil Null behavior α, X.B Bind ω(α, X).B Weak split ρ(L).B Strong split ǫ(t).B Delay B + B Choice (B)i Iteration Θ Thanatos δ(t, B) Duration K Process name P ∈ 3DP - 3D process S[B] S ∈ S, B ∈ B Pa, XP Composed 3D process N ∈ N - 3D network nil Empty network P P ∈ 3DP NN Parallel 3D processes

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Shape Calculus

S0[B0] S1[B1] S0[B′

0]a, Y S1[B′ 1]

a, X1 a, X0 b, X ′

1

c, X ′ S0[B′′

0 ]

S1[B′′

1 ]

a, Y

P0, P1, Y t ρ ∆

Evolution of 3D processes in the Shape Calculus. Processes S0[B0] and S1[B1] are involved in a bind on the channel a, Y and in a subsequent strong split.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Shape Calculus model for BR

• Cells communicate

each other via direct contact

• RANK/RANKL/OPG

are surface-bound proteins

• Cells as 3D Processes,

communicating by binding

• Biochemical signals as

channels exposed at the surface of the 3D process

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Shape Calculus model for BR

• Cells communicate each
• ther via direct contact
• RANK/RANKL/OPG are

surface-bound proteins

• Cells as 3D Processes,

communicating by binding

• Biochemical signals as channels

exposed at the surface of the 3D process

RANKL/OPG signaling. (pre-osteoblast, mature osteoblast) → (mature osteoblast, active osteoblast)

ρ(rankl, X1) rankl, X2 rankl, Y ρ(rankl, X2) BOb rankl, X2 rankl, X1 SOb SOb SOb SOb SOb SOb

SOb[rankl, X1.ρ(rankl, X1).BOPG] SOb[rankl, X2.ρ(rankl, X2).BOb] ↓ SOb[ρ(rankl, X1).BOPG]rankl, YSOb[ρ(rankl, X2).BOb] ↓ SOb[BOPG] SOb[BOb]

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

BR specification

Tissue, BMU Tissue

def

= (ABMUi)a

i=1 (QBMUj)q j=1

ABMU

def

= (Oyi)

nOy i=1 (Ocj)nOc j=1 (Obk)nOb k=1

QBMU

def

= (Oyi)

nOy i=1

Osteocyte Oy

def

= SOy[(can, X + can, X)kOy + consume, X.Θ] Osteoclast Oc

def

= SOc[BOc] BOc

def

= δ(tOc, (consume, X.resorb)∞).Θ.mineral, XkOc Osteoblast Ob

def

= SOb[rank, X.ρ(rank, X).BOPG + ǫ(tPb).BOPG] BOPG

def

= rank, X.ρ(rank, X).BOb + BOb BOb

def

= δ(tOb, (mineral, X.form)∞).Θ

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Encoding of 3D processes as agents

• Stochastic actions Binding affected by the rates of the

actions involved (rated output/passive input)

• Perception distance Agents can communicate within a

given radius

• Prototype in Repast Simphony: support for Shape

Calculus models (simplified shapes), stochastic simulation implemented through discrete-event scheduler

A B a, w2 a, w1 a, λ

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Encoding of 3D processes as agents

• Stochastic actions Binding affected by the rates of the

actions involved (rated output/passive input)

• Perception distance Agents can communicate within a

given radius

• Prototype in Repast Simphony: support for Shape

Calculus models (simplified shapes), stochastic simulation implemented through discrete-event scheduler

A

(a,λ)d

− − − → B

a,w1

− − → D

a,w3

− − → C

a,w2

− − → d

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Encoding of 3D processes as agents

• Stochastic actions Binding affected by the rates of the

actions involved (rated output/passive input)

• Perception distance Agents can communicate within a

given radius

• Prototype in Repast Simphony: support for Shape

Calculus models (simplified shapes), stochastic simulation implemented through discrete-event scheduler

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

BR model

Motion

• Output channels release “molecule concentrations”

(RANKL and Oc death factors) diffusing according to a CA-like rule

• Cells move according to a biased random walk, influenced

by compatible molecular bias

Parametrization

• cells number, lifetime, size, resorption and formation rate

taken from experimental works

• diffusion coefficients and perception radii tuned for having

realistic remodelling times

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Results

Experimental evidence: RANKL concentration inversely related to bone turnover and density. Aging affects bone structure.

Two parameter configurations

• healthy, with regular RANKL production and cellular activity;
• pathological, with an overproduction of RANKL and a reduced cellular

activity (i.e. aging).

Healthy Osteoporotic Stdev H Stdev O

time [days] BMD [mg/cm2]

40 simulations of agent-based model, 1 remodelling cycle

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Results

Experimental evidence: RANKL concentration inversely related to bone turnover and density. Aging affects bone structure.

Two parameter configurations

• healthy, with regular RANKL production and cellular activity;
• pathological, with an overproduction of RANKL and a reduced cellular

activity (i.e. aging). Normal Osteopenia Osteoporosis After 7 remodelling cycles (∼ 7 years), pathological patient goes osteoporosis.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Results

H O H O H O

1st cycle 2nd cycle 3rd cycle

Positioning of signalling osteocytes affect bone microstructure

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Conclusion

• Extension of Shape Calculus
• New agent-based, stochastic and spatial model of BR, based
• n formal language specification, and able to reproduce

normal and defective dynamics

• Too many parameters if using the full expressive power of the
• calculus. Practically, simplifications are needed.

[1]

• N. Paoletti, P. Li`
• , E. Merelli and M. Viceconti.

Multi-level Computational Modeling and Quantitative Analysis of Bone Remodeling. In IEEE/ACM TCBB, 9(5), pp. 1366-1378, 2012. [2]

• N. Paoletti, P. Li`
• , E. Merelli and M. Viceconti.

Osteoporosis: a multiscale modeling viewpoint. In CMSB 2011, 9(5), ACM, pp. 183-193, 2011. [3]

• P. Li`
• , E. Merelli, N. Paoletti and M. Viceconti.

A combined process algebraic and stochastic approach to bone remodeling. In CS2Bio 2011, ENTCS 277, pp. 41-52, 2011.

1.2- Formal Analysis of Bone Pathologies

Hybrid Ap- proximation CTMC ODE Probabilistic Model Checking Stabili- zation Simulation Sensitivity Analysis Model Checking Parameter Synthesis Osteopo- rosis Osteomy- elitis

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

ODE model

From [Komarova et al. Bone 33.2 (2003): 206-215] (Osteoclasts) ˙ x1 =α1x1

g11x2 g21 − β1x1

(Osteoblasts) ˙ x2 =α2x1

g12x2 g22 − β2x2

(Bone density) ˙ z = − k1x1 + k2x2

5 10 15 20 5 10 Time (d) Osteoclasts 200 400 50 100 Time (d) Osteoblasts 200 400 −40 −20 Time (d) ”Bone”

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Calibration

Aim: Tune ODE parameters (growth and death rates) so that max

• num. of bone cells agree with more recent evidence (max 10 Oc

and 100 Ob, against 20 Oc and 2000 Ob of the original model).

Workflow

• Local sensitivity of bone cells
• Global sensitivity of max Oc and Ob
• Fitting and estimation against (made up) observations

10 20 −6 −4 −2 2 Time (d) Oc sensitivity α1 β1 α2 β2 10 20 −100 100 Time (d) Ob sensitivity α1 β1 α2 β2

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Calibration

Aim: Tune ODE parameters (growth and death rates) so that max

• num. of bone cells agree with more recent evidence (max 10 Oc

and 100 Ob, against 20 Oc and 2000 Ob of the original model).

Workflow

• Local sensitivity of bone cells
• Global sensitivity of max Oc and Ob
• Fitting and estimation against (made up) observations
• ●
• ●
• ●
• ●
• 0.0

0.2 0.4 0.6 0.8 1.0 50 100 150 200

Osteoblasts

alpha1 Maximum cell number

• ●
• ●
• ●
• 0.5

0.6 0.7 0.8 0.9 1.0 50 100 150 200

Osteoblasts

beta1 Maximum cell number

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Calibration

Aim: Tune ODE parameters (growth and death rates) so that max

• num. of bone cells agree with more recent evidence (max 10 Oc

and 100 Ob, against 20 Oc and 2000 Ob of the original model).

Workflow

• Local sensitivity of bone cells
• Global sensitivity of max Oc and Ob
• Fitting and estimation against (made up) observations

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

CTMC model

Derived from ODE model and implemented in the model checker PRISM

Osteoclasts: [ ] 0 < x1 < xmax

1

∧ x2 > 0 → α1xg11

1

xg21

2

: x1 = x1 + 1 [ ] x1 > 0 → β1x1 : x1 = x1 − 1 [resorb] x1 > 0 → k1x1 : true Osteoblasts: [ ] 0 < x2 < xmax

2

∧ x1 > 0 → α2xg12

1

xg22

2

: x2 = x2 + 1 [ ] x2 > 0 → β2x2 : x2 = x2 − 1 [form] x2 > 0 → k2x2 : true Bone reward: [resorb] true : 1 [form] true : 1

100 200 300 400 20 40 60 80 100

Osteoblasts

Time [days] Number of cells 5 10 15 20 2 4 6 8 10

Osteoclasts

Time [days] Number of cells

Bone

Relative Density Time [days]

0 100 200 300 20

• 20
• 40
• 60
• 80

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Hybrid approximation

• Optimal approximation of the non-linear ODEs into a

piecewise multiaffine (PMA) system

• Overapproximation into a transition system where states are

hyper-rectangles (reachable sets) and transitions are possible trajectories between reachable sets. Analyzed with the tool RoVerGeNe by Batt et al.

˙ x1 = α1

ns1+1

• i=1

r(x1, θ(1)

i

, θ(1)

i+1, y(1) i

, y(1)

i+1)) ns2+1

• i=1

r(x2, θ(2)

i

, θ(2)

i+1, y(2) i

, y(2)

i+1)) − β1x1

• r is the ramp function

r(x, θ1, θ2, y1, y2) =      y2 (x ≥ θ2) y1 + (y2 − y1) (x−θ1)

(θ2−θ1)

(θ1 ≤ x < θ2) y1 (x < θ1)

• θ(k)

1

= xmin

k

< θ(k)

2

< · · · < θ(k)

ns1+1 = xmax k

are the segments identified by the approximation algorithm for xk

• y(k)

i

is the ith sampled point for xk

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Hybrid approximation

• Optimal approximation of the non-linear ODEs into a

piecewise multiaffine (PMA) system

• Overapproximation into a transition system where states are

hyper-rectangles (reachable sets) and transitions are possible trajectories between reachable sets. Analyzed with the tool RoVerGeNe by Batt et al.

By the convexity property of multi-affine funs, we just look at the vertices

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Stabilization properties

• Stabilisation is usually the existence of a unique fixpoint

(state) or attractor (region) that is always eventually reached: in LTL, FG(prop).

• In our case, we follow the pattern φ =

⇒ G(prop), focusing

• n particular regions (precondition φ) of the state space
• Note: with conservative extension of PMA, we can be sure
• nly of positive answers (simulation weakly preserves LTL

SAT relation)

Down-regulation of Oc by Ob

With high numbers of Obs, Ocs never grow

• PMA model √

(x2 > θ(2)

3 ) → G(∧11 i=2(x1 < θ(1) i

∧x2 > θ(2)

3 ) =

⇒ X(x1 < θ(1)

i

))

• CTMC model √, ∀θ(1) = 1, . . . , xmax

1

P=1 [G((x1 < θ(1) ∧ x2 > θ(2)

3 ) =

⇒ X(x1 < θ(1)))]

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Stabilization properties

• Stabilisation is usually the existence of a unique fixpoint

(state) or attractor (region) that is always eventually reached: in LTL, FG(prop).

• In our case, we follow the pattern φ =

⇒ G(prop), focusing

• n particular regions (precondition φ) of the state space
• Note: with conservative extension of PMA, we can be sure
• nly of positive answers (simulation weakly preserves LTL

SAT relation)

Converse

With low numbers of Obs, Ocs never decay

• PMA model √

(x2 < θ(2)

2 ) =

⇒ G(∧11

i=4(x1 > θ(1) i

∧x2 < θ(2)

2 ) =

⇒ X(x1 > θ(1)

i

))

• CTMC model √, ∀θ(1) = 1, . . . , xmax

1

P=1 [G((x1 > θ(1) ∧ x2 < θ(2)

1 ) =

⇒ X(x1 > θ(1)))]

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Parameter synthesis

• We find the regions of the parameter space s.t. Ocs and Obs

have a fixed upper bound ((x1 < θ(1)

10 ) ∧ (x2 < θ(2) 20 )) → G((x1 < θ(1) 10 ) ∧ (x2 < θ(2) 20 ))

• Note: synthesis = estimation/identification (where a

single param value is provided)

• Synthesis method (Batt et al.), based on hierarchical

decomposition of param space guided by affine inequalities

• ver the parameters (termination ensured by the existence of

param equivalence classes).

3.73 4.54 4.65 8 ·10−2 1.6 2.04 2.4 2.8 β2 β1

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Robust stabilisation wrt init conditions

Stabilization (bone homeostasis) reached regardless init num of

• steoclasts

ODE model, global sensitivity analysis

CTMC model, probabilistic (quantitative) model checking, parametric analysis

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Defective Dynamics

Osteoporosis: characterized by low bone mass and structural fragility, results from ageing and defects in the RANK/RANKL/OPG pathway. Osteomyelitis: a bacterial infection of bone (S. aureus) that affects the signalling level, and rapidly leads to bone loss and necrosis.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Defective Dynamics - Extended ODEs

(Osteoclasts) ˙ x1 =α1x

g11(1+f11 B

s )

1

x

g21(1−f21 B

s )+gpor

2

− kageingβ1x1 (Osteoblasts) ˙ x2 =α2x

g12/(1+f12 B

s )

1

x

g22−f22 B

s

2

− kageingβ2x2 (S. aureus) ˙ B =(γB − V )B · log( s B )

• B follows a logistic growth
• fij: effects of the infection on auto- and para- regulation

factors gij

• V : medical treatment (V > γB, bacteriocide; V = γB,

bacteriostatic)

• gpor: effectiveness of RANKL expression
• kageing: ageing factor increasing cells death rate

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Defective Dynamics - Extended ODEs

(Osteoclasts) ˙ x1 =α1x

g11(1+f11 B

s )

1

x

g21(1−f21 B

s )+gpor

2

− kageingβ1x1 (Osteoblasts) ˙ x2 =α2x

g12/(1+f12 B

s )

1

x

g22−f22 B

s

2

− kageingβ2x2 (S. aureus) ˙ B =(γB − V )B · log( s B )

BACTERIOSTATIC TREATMENT BACTERIOCIDE TREATMENT

Simulation of bacteriostatic (V = γB) and antibiotic (V > γB) treatments.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Probabilistic model checking results







  

Stochastic vs ODE model

a) b) c)

Model checking-based clinical estimator: density change rate.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Conclusion

• Formal analysis gives not only biological insights

(stabilization), but also estimators of clinical importance (bone diseases)

• Stochastic model closer to original ODEs → Probabilistic

model checking results are “to be trusted”

• Unlike hybrid approximation, no effective methods for

parameter synthesis of CTMCs

[1]

• E. Bartocci, P. Li`
• , E. Merelli and N. Paoletti.

Multiple verification in complex biological systems: the bone remodelling case study. In TCSB XIV, LNCS 7625, pp. 53-76, 2012. [2]

• P. Li`
• , N. Paoletti, M.A. Moni, K. Atwell, E. Merelli and M. Viceconti.

Modelling osteomyelitis. In BMC Bioinformatics, 13(Suppl 14), 2012. [3]

• P. Li`
• , E. Merelli and N. Paoletti.

Multiple verification in computational modeling of bone pathologies. In CompMod2011, EPTCS 67, pp. 82-96, 2011.

2- Hybrid Approach to the Scheduling of Multiple Therapies

D − CGF (extension

• f CGF PA)

Controlled Switched Systems Stochastic Hybrid Automata Optimal Control STI (Ther- apy Schedul- ing)

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Optimal Scheduling of Multiple Therapies

Aim

Structured Therapy Interruption (STI) is the programmed interruption of a medication for a period of time. STI Problem. Find the “best” combination of multiple on-off therapies to recover a dynamic disease model CGF − → ODE systems

• Based on CGF (Chemical Ground Form) (Cardelli, TCS

2008, pp. 261–281), a subset of stochastic CCS

• (Bortolussi et al., Mathematics in Computer Science 2(3), pp.

465–491, 2009): hybrid automata semantics for CGF

• D-CGF (Disease CGF) extends CGF to describe complex

disease processes, with a hybrid dynamical systems semantics.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Optimal Scheduling of Multiple Therapies

Aim

Structured Therapy Interruption (STI) is the programmed interruption of a medication for a period of time. STI Problem. Find the “best” combination of multiple on-off therapies to recover a dynamic disease model CGF − → Hybrid automata

• Based on CGF (Chemical Ground Form) (Cardelli, TCS

2008, pp. 261–281), a subset of stochastic CCS

• (Bortolussi et al., Mathematics in Computer Science 2(3), pp.

465–491, 2009): hybrid automata semantics for CGF

• D-CGF (Disease CGF) extends CGF to describe complex

disease processes, with a hybrid dynamical systems semantics.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Optimal Scheduling of Multiple Therapies

Aim

Structured Therapy Interruption (STI) is the programmed interruption of a medication for a period of time. STI Problem. Find the “best” combination of multiple on-off therapies to recover a dynamic disease model D-CGF − → Hybrid dynamical systems

• Based on CGF (Chemical Ground Form) (Cardelli, TCS

2008, pp. 261–281), a subset of stochastic CCS

• (Bortolussi et al., Mathematics in Computer Science 2(3), pp.

465–491, 2009): hybrid automata semantics for CGF

• D-CGF (Disease CGF) extends CGF to describe complex

disease processes, with a hybrid dynamical systems semantics.

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Syntax

Two kind of processes:

• Individuals: standard CGF processes that represent the population

in the disease model (e.g. susceptible, infected, . . . ). Interpreted as continuous variables in the hybrid semantics.

• Therapies: processes for modelling discrete interventions on the

disease scenario (drug dosage). Interpreted as discrete switches in the semantics.

E ::= 0 | X = M, E Reagents M ::= 0 | π.P + M Molecule P ::= 0 |(XP) Solution π ::= τ r | ?xr | !xr Actions (r ∈ R+) D ::= 0 | Y = C, D Set of therapy definitions C ::= 0 | π.T + C Therapy T ::= 0 |(Y T) Combination of therapies D − CGF ::= (E, P, D, T) A D-CGF model

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Syntax

Two kind of processes:

• Individuals: standard CGF processes that represent the population

in the disease model (e.g. susceptible, infected, . . . ). Interpreted as continuous variables in the hybrid semantics.

• Therapies: processes for modelling discrete interventions on the

disease scenario (drug dosage). Interpreted as discrete switches in the semantics.

E ::= 0 | X = M, E Reagents M ::= 0 | π.P + M Molecule P ::= 0 |(XP) Solution π ::= τ r | ?xr | !xr Actions (r ∈ R+) D ::= 0 | Y = C, D Set of therapy definitions C ::= 0 | π.T + C Therapy T ::= 0 |(Y T) Combination of therapies D − CGF ::= (E, P, D, T) A D-CGF model

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Epidemics with therapies

We consider a SIR epidemic model + therapies:

• D1, a vaccination (makes immune the susceptible

population)

• D2, an antiviral (increases the recovery rate)

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

SIR + therapies D-CGF model

We build on a open population SIR (births and deaths):

D-CGF model M = (E, P, D, T)

E: S = τ b

S1.(SS) + τ µ S2.0 + ?jρ.R + ?iβ.I

I = τ b

I1.(IS) + τ µ I2 .0 + ?hk.R + !iβ.I + τ ν I3.R

R = τ b

R1.(RS) + τ µ R2.0

D: D1off = τ r1on

1on .D1on

D1on = τ r1off

1off .D1off + !jρ.D1on

D2off = τ r2on

2on .D2on

D2on = τ r2off

2off .D2off + !hk.D2on

C : D1off D2off

List of reactions

S + D1on →ρ R + D1on Susceptible becomes recovered with therapy 1 I + D2on →k R + D2on Infected becomes recovered with therapy 2 D1off →r1on D1on D1 switched on D1on →r1off D1off D1 switched off D2off →r2on D2on D2 switched on D2on →r2off D2off D2 switched off

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

SIR + therapies D-CGF model

We build on a open population SIR (births and deaths):

D-CGF model M = (E, P, D, T)

E: S = τ b

S1.(SS) + τ µ S2.0 + ?jρ.R + ?iβ.I

I = τ b

I1.(IS) + τ µ I2 .0 + ?hk.R + !iβ.I + τ ν I3.R

R = τ b

R1.(RS) + τ µ R2.0

D: D1off = τ r1on

1on .D1on

D1on = τ r1off

1off .D1off + !jρ.D1on

D2off = τ r2on

2on .D2on

D2on = τ r2off

2off .D2off + !hk.D2on

C : D1off D2off

List of reactions

S + D1on →ρ R + D1on Susceptible becomes recovered with therapy 1 I + D2on →k R + D2on Infected becomes recovered with therapy 2 D1off →r1on D1on D1 switched on D1on →r1off D1off D1 switched off D2off →r2on D2on D2 switched on D2on →r2off D2off D2 switched off

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

SIR + therapies D-CGF model

We build on a open population SIR (births and deaths):

D-CGF model M = (E, P, D, T)

E: S = τ b

S1.(SS) + τ µ S2.0 + ?jρ.R + ?iβ.I

I = τ b

I1.(IS) + τ µ I2 .0 + ?hk.R + !iβ.I + τ ν I3.R

R = τ b

R1.(RS) + τ µ R2.0

D: D1off = τ r1on

1on .D1on

D1on = τ r1off

1off .D1off + !jρ.D1on

D2off = τ r2on

2on .D2on

D2on = τ r2off

2off .D2off + !hk.D2on

C : D1off D2off

List of reactions

S + D1on →ρ R + D1on Susceptible becomes recovered with therapy 1 I + D2on →k R + D2on Infected becomes recovered with therapy 2 D1off →r1on D1on D1 switched on D1on →r1off D1off D1 switched off D2off →r2on D2on D2 switched on D2on →r2off D2off D2 switched off

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Derivation of hybrid dynamical systems

We focus on Controlled Switched Systems (CSS), a class of hybrid dynamical systems where the external control input is given by the discrete operation mode.

System ˙ x = f (x, q) y = g(x, q) Controller

Drug scheduler

y

Observables (e.g. pa- tient medical data)

q

Therapy (combi- nation of drugs)

Control loop (x continuous state, y observable output, q discrete mode)

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Derivation of hybrid dynamical system

For a correct derivation of the semantics, the D-CGF model must be well-formed, i.e. informally, its drug terms can be interpreted as discrete variables.

1 Build the SD-graph, a digraph whose nodes are therapy

terms in D and arcs connect couples (Yi, Yj) involved in a switch, i.e. an action π that consumes Yi, and produces Yj: ∆(π, Yj) = −1 and ∆(π, Yj) = 1

2 Consider its connected components

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Derivation of hybrid dynamical system

For a correct derivation of the semantics, the D-CGF model must be well-formed, i.e. informally, its drug terms can be interpreted as discrete variables.

1 Build the SD-graph, a digraph whose nodes are therapy

terms in D and arcs connect couples (Yi, Yj) involved in a switch, i.e. an action π that consumes Yi, and produces Yj: ∆(π, Yj) = −1 and ∆(π, Yj) = 1

2 Consider its connected components

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Derivation of hybrid dynamical system

3 Derivation of the discrete modes from the SD-graph 4 Modified rate vector:

For each mode qi, φqi is obtained according to the drug terms active in qi:

. . . τ1on τ1off τ2on τ2off j h φ = . . . T1offr1on T1onr1off T2offr2on T2onr2off T1onρS T2onkI φq = q1 . . . r1on r2on q2 . . . r1off r2on ρS q3 . . . r1on r2off kI q4 . . . r1off r2off ρS kI

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Derivation of hybrid dynamical system

5 Controlled switched system given, at each mode qi, by the

product of the stochiometric matrix (restricted to i) and the rate vector at qi ˙ x(qi) = M|E · φqi

˙ x(q1) =    ˙ S = b(S + I + R) − βSI − µS ˙ I = βSI − µI − νI ˙ R = νI − µR ˙ x(q2) =    ˙ S = b(S + I + R) − βSI − µS − ρS ˙ I = βSI − µI − νI ˙ R = νI − µR + ρS ˙ x(q3) =    ˙ S = b(S + I + R) − βSI − µS ˙ I = βSI − µI − (ν + k)I ˙ R = (ν + k)I − µR ˙ x(q4) =    ˙ S = b(S + I + R) − βSI − µS − ρS ˙ I = βSI − µI − (ν + k)I ˙ R = (ν + k)I − µR + ρS

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Non-linear optimal control definition

We define the following Model Predictive Control (MPC) problem to compute the optimal administration of therapies.

min

q t+Tp

• t

{Rq(k)1 + Qx(k)1} dk

• subj. to

x(t) ∈ [0, 1]3 q(t) ∈ {0, 1}4

• i qi(t) = 1

˙ x(t) = f (x(t), q(t)) x(0) = x0 x(Tf ) ∈ Xf ∀t ∈ [0, Tf ]

• Tp is the prediction horizon
• ∆t = 1/365 is the discrete time

step (1 day)

• q: controlled modes; x:

continuous state

• Xf = {[S I R]T ∈

[0, 1] × [0, I0] × [0, 1]} is set of terminal states (the final num. of infected will not be higher than the init one)

• R and Q are the

weights/penalties on the controlled inputs and on the system state

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Controlling non-linear hybrid systems

No exact algorithms for the optimal control of generic non-linear hybrid systems (except for specific classes, up to PMA)

Embedding approach

• Original system (SOCP) is embedded into a larger class of

systems without discontinuities (EOCP)

• Relaxes mode vector q ∈ {0, 1}n into ˜

q ∈ [0, 1]n

• The optimal continuous modes of the EOCP are projected

back as discrete modes

• Facts:
• SOCP has solution =

⇒ EOCP has solution

• EOCP solution can be approximated with arbitrary precision

to SOCP solution

• Mode projection is theoretically justified and practically

effective to approximate EOCP solution

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Embedded problem

min

˜ q t+Tp

• t

{˜ qi(k)(R·,i1 + Q˜ x(k)1)} dk

• subj. to

˜ x(t) ∈ [0, 1]3 ˜ q(t) ∈ [0, 1]4

• i ˜

qi = 1 ˙ ˜ x(t) =

i ˜

qi · fqi (˜ x(t)) ˜ x(0) = x0 ˜ x(Tf ) ∈ Xf ∀t ∈ [0, Tf ]

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Application to the H1N1 Case Study

We start from a SIR model of H1N1 influenza occurred in US in 2009 (Feng et al. AAPS, 13(3):427–437, 2011)

• Large R −

→ drug dosage tends to be avoided (severe side effects or early stage of the disease)

• Small R −

→ more prominent drug dosage (little side effects or mature stage of disease)

Treatment scenarios

1 Low penalties to D1 and D2 2 High penalty to D1 (antiviral favoured) 3 High penalty to D2 (vaccination favoured) 4 High penalty to the combination of D1 and D2

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Results

Evolution of epidemics and administration of vaccination and antiviral treatments

1 Drugs

Scenario 1

D1 D2 0.5 1 Population 1 Drugs

Scenario 2

50 100 150 200 250 300 350 0.5 1 Time (d) Population

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Results

100 200 300 20 40 60 80

4.85 86.5 4.86 4.87

Time (d) Cumulative Infections Scenario 1 Scenario 2 Scenario 3 Scenario 4

Cumulated cases of infection under the four treatment scenarios

100 200 300 100 200 300 400

52.66 367.46 55.01 54.85

Time (d) Cumulated cost Scenario 1 Scenario 2 Scenario 3 Scenario 4

Cumulated costs of optimal strategies

Formal Computational Modelling of Bone Physiology and Disease Processes

• N. Paoletti

Conclusions

• From (stochastic) process algebra to hybrid non-linear
• ptimal control
• Approximated methods effective for the control of non-linear

hybrid systems

• Directions: game-theoretic extension, stochasticity in

controlled system, . . . [1]

• P. Li`
• , E. Merelli and N. Paoletti.

Disease processes as hybrid dynamical systems. In of HSB 2012, EPTCS 92, pp. 152-166, 2012.