Bayesian inference in dynamic modeling of biological systems (a - PowerPoint PPT Presentation

Bayesian inference in dynamic modeling of biological systems (a modeling practitioner’s perspective) Karl Thomaseth National Research Council Institute of Biomedical Engineering, ISIB-CNR Padova, Italy

3 H-Glucose - 3 H-H 2 O Kinetics

Early Estimation Approach • Nonlinear least squares n � ˆ w i ( y i − f ( t i , ϑ )) 2 ϑ = arg min ϑ i =1 • Justification y i = f ( t i , ϑ ) + e ( t i ); e ( t i ) ∼ N (0 , σ 2 ( t i )) n ( y i − f ( t i , ϑ )) 2 log p ( y | ϑ ) = const. − 1 � σ 2 ( t i ) 2 i =1 ⇒ ˆ ϑ = ˆ if w i = σ − 2 ( t i ) (known !) = ϑ ( ML ) .

Computations • Restricted step Gauss-Newton, aka Levenberg-Marquardt n � ∇ ϑ log p ( y | ϑ ) = w i ( y i − f ( t i , ϑ )) ∇ ϑ f ( t i , ϑ ) i =1 n � ∇ 2 w i ∇ ϑ f ( t i , ϑ ) ∇ ϑ f ( t i , ϑ ) ′ ϑ 2 log p ( y | ϑ ) = − (+ . . . ) i =1 • Iterations ϑ k +1 = ˆ ˆ ϑ k − [ ˜ ∇ 2 ϑ 2 log p ( y | ϑ k ) + λI ] − 1 ∇ ϑ log p ( y | ϑ k ) • Estimator variance (Cramer-Rao LB) � − 1 ⇒ − [ ˜ � var (ˆ ϑ 2 log p ( y | ˆ ϑ )] − 1 E [ −∇ 2 ∇ 2 ϑ ) ≥ ϑ 2 log p ( y | ϑ )]

Maximum a Posteriori • Prior - Posterior distributions ∼ p ( ϑ ) N ( ϑ 0 , Σ 0 ) p ( ϑ | y ) ∝ p ( y | ϑ ) p ( ϑ ) � n � ( y i − f ( t i , ϑ )) 2 − 1 � + ( ϑ − ϑ 0 ) ′ Σ − 1 log p ( ϑ | y ) 0 ( ϑ − ϑ 0 ) = + c. σ 2 ( t i ) 2 i =1 • Implications n � w i ( y i − f ( t i , ϑ )) ∇ ϑ f ( t i , ϑ ) − Σ − 1 ∇ ϑ log p ( ϑ | y ) 0 ( ϑ − ϑ 0 ) = i =1 � n � � w i ∇ ϑ f ( t i , ϑ ) ∇ ϑ f ( t i , ϑ ) ′ + Σ − 1 ˜ ∇ 2 ϑ 2 log p ( ϑ | y ) − = < 0 0 i =1

Dynamic Modeling Classes • Regression (black box) models, e.g. f ( t, ϑ ) = Ae − α t + Be − β t • Comprehensive models based on physi(ologi)cal principles – conservation of mass ⇒ compartmental models – conservation of energy (chemical, mechanical, thermal) ⇒ multiphysics models – anatomy (circulation), metabolic pathways • Minimal models – highly aggregated physiological models – identifiability of parameters

PBPK - Models

Body Fluid Compartments 7.5% Dense CT - Cartilage 20% 7.5% 55% 2.5% Interstitial- Plasma Intracellular Transcellular Lymph 7.5% Bone

Two Compartment Approximation 7.5% Dense CT - Cartilage Extracellular Intracellular 20% 7.5% 55% 2.5% Interstitial- Plasma Intracellular Transcellular Lymph 7.5% Bone

Model Equivalence A B u y u y k12 Cl12 1 2 1 2 k21 k01 k02 Cl1 Cl2 C D Qc Qc Cart Cart Q1 1 Q1 C1 1 C1 Cl1 Cl1 Q2 2 Q2 C2 2 Cl2 C2 Q3 kidney Cl2 C3 Cl3

Bayesian inference: Example Modeling Population Kinetics of Free Fatty Acids in Isolated Rat Hepatocytes using Markov Chain Monte Carlo Alessandra Pavan, Karl Thomaseth Institute of Biomedical Engineering, Padova, Italy Anna Valerio Department of Clinical and Experimental Medicine, University of Padova, Italy

Cellular Fuels - Glucose BLOOD Glycogen Adipocyte Glucose Glucose 6-P Glucose Liver – Acetyl CoA TCA Pyruvate Acetyl CoA Ketone Ketone bodies Triacylglycerol TCA bodies HSL Fatty acids Glycerol Fatty acids – ATP Fatty acids Fatty acids Gluconeogenic Glycerol precursors Ketone Amino acids bodies Fatty acids Insulin Ketone bodies + Amino Protein + Pancreas acids Glucose Glucose Acetyl CoA Muscle TCA

Cellular Fuels - Free Fatty Acids BLOOD Glycogen Adipocyte Glucose 6-P Glucose Glucose Liver – Acetyl CoA TCA Pyruvate Acetyl CoA Ketone Ketone bodies Triacylglycerol TCA bodies HSL Fatty acids Glycerol Fatty acids – ATP Fatty acids Fatty acids Gluconeogenic Glycerol precursors Ketone Amino acids bodies Fatty acids Insulin Ketone bodies + Amino Protein + Pancreas acids Glucose Glucose Acetyl CoA Muscle TCA

Cellular Fuels - Control BLOOD Glycogen Adipocyte Glucose Glucose 6-P Glucose Liver – Acetyl CoA TCA Pyruvate Acetyl CoA Ketone Ketone bodies Triacylglycerol TCA bodies HSL Fatty acids Glycerol Fatty acids – ATP Fatty acids Fatty acids Gluconeogenic Glycerol precursors Ketone Amino acids bodies Fatty acids Insulin Ketone bodies + Amino Protein + Pancreas acids Glucose Glucose Acetyl CoA Muscle TCA

Hormonal Control of Ketogenesis adipose tissue blood liver FFA FFA triglyceride ketone bodies triglycerides phospholipids glucagon & epinephrine insulin

Experimental Protocol N 9 Oleate Dose 8 0 mmol/L (N=30) 0.25 mmol/L (N=18) 7 0.5 mmol/L (N=23) 6 1.0 mmol/L (N=13) 5 4 3 2 1 0 1 2 3 4 5 6 Protocol Oleate (FFA) no yes yes no yes yes Epinephrine no no yes yes no yes Insulin no no no no yes yes

Ketone Body Production FFA �ß - Oxidation k 2 k 1 k 3 k 4 AcAc Dose Dose k 3 ≈ k 2 k 2 FFA AcAc BOH FFA + k 4 BOH k 1 k 1

Experimental Data - Examples D=0 D=0.25 D=0.5 D=1 1 0.4 1.5 2 mmol/L 1 0.2 0.5 1 0.5 0 0 0 0 0 50 100 0 50 100 0 50 100 0 50 100 1 0.4 1.5 2 mmol/L 1 0.5 0.2 1 0.5 0 0 0 0 0 50 100 0 50 100 0 50 100 0 50 100 0.4 1 2 1 mmol/L 0.2 0.5 0.5 1 0 0 0 0 0 50 100 0 50 100 0 50 100 0 50 100 2 0.2 0.4 1 mmol/L 1 0.1 0.2 0.5 0 0 0 0 0 50 100 0 50 100 0 50 100 0 50 100 min min min min

Minimal FFA – KB Kinetic Model • System dynamics with known + random FFA dose  − ( k 1 + k 2 ) x 1 ( t ) ; x 1 ( t ) ˙ = x 1 (0) = D + b  x 2 ( t ) ˙ = 4 k 2 x 1 ( t ) ; x 2 (0) = 0   y 1 ( t ) = x 1 ( t j ) + x 1 b + ǫ 1 ( t ) ; t = { t j } 1 ,...,n  y 2 ( t ) = x 2 ( t ) + x 2 b + p b t + ǫ 2 ( t )  • Constraints and parameterization k = k 1 + k 2 ( k > 0) ⇒ ϑ k = log k (0 ≤ ϕ ≤ 1) ⇒ ϑ ϕ = log( ϕ/ (1 − ϕ )) ϕ = k 2 / ( k 1 + k 2 ) ⇒ ϑ = [log k, logit ϕ, log p b , log x 1 b , log x 2 b ] • Effects of covaritates ϑ [ k,ϕ,p b ,x 1 b ,x 2 b ] = X · µ [ k,ϕ,p b ,x 1 b ,x 2 b ] ; X = [1 , ∆ N cells , D, I Ins , I Epi ]

Hierarchical Population Kinetic Model 3. Prior distributions Hyperparameters p ( µ, Σ , σ 2 ) = p ( µ ) p (Σ) p ( σ 2 ) 2 µ Σ σ 2. Between-individual variability ϑ i ∼ p ( µ, X i , Σ) X i θ i ∼ N ( X i · µ, Σ) 1. Within-individual variability p ( y | ϑ i , D i , t ij , σ 2 ) ∼ y ij y ij t ij D i N ( f ( ϑ i , D i , t ij ) , σ 2 ) ∼

Bayesian Inference • Joint Posterior Probability �� n i � � N π ( ϑ, µ, σ 2 , Σ j p ( y ij | ϑ i , σ 2 ) | y ) ∝ p ( ϑ i | µ, Σ) × i � �� � x p ( µ ) p ( σ 2 ) p (Σ) × • Paradigm: inference, e.g. mean, SD, median, 95% CI, based on marginal distributions: π ( ϑ | y ) , π ( µ | y ) , . . . � g ( x ) π ( x | y ) d x = E π ( g ( x )) χ • Monte Carlo approach: generate iid random samples x 1 , ..., x n , x i ∼ π ( x | y ) ∀ i = 1 , . . . , n n � 1 g ( x i ) − → E π ( g ( x )) n i =1

Markov Chain Monte Carlo • Problem: the joint posterior distribution is impossible to obtain analytically in general (exception conjugate priors) • Solution: generate, non-independent, Markov chains of random variables sampled from the full conditional distributions using: – Gibbs sampler – Metropolis Hastings – Mixed algorithms: Gibbs + Metropolis • Computationally intensive + need to discard initial samples (burnin) to reach steady-state distributions

Gibbs Sampler: x ∼ π ( x | y ) x t = ( x 1 , x 2 , · · · , x k ) ↓ x 1 ∼ π ( x 1 | x 2 , x 3 , · · · x k , y ) ⇓ x 2 ∼ π ( x 2 | x 1 , x 3 , · · · x k , y ) ⇓ x 3 ∼ π ( x 3 | x 1 , x 2 , · · · x k , y ) . . . x k ∼ π ( x k | x 1 , x 2 , · · · x k − 1 , y ) ↓ x t +1 = ( x 1 , x 2 , · · · , x k ) ∼ π ( x 1 , x 2 , · · · x k − 1 , x k | y ) & x i ∼ π ( x i | y ) i = 1 , . . . , k

Metropolis Hastings: x i ∼ π ( x i | x ( − i ) , y ) • Applicable if full conditional distributions are not standard – given x t i – generate candidate value w i ∼ q ( ·| x i ) – accept ( x t +1 = w i ) with probability i � � 1 , π ( w i | ... ) q ( x i | w i ) α ( w i | x i ) = min π ( x i | ... ) q ( w i | x i ) • Random Walk Metropolis ( q ( x i | w i ) = q ( w i | x i ) ) – w i = x t i + δ i δ i ∼ N (0 , Σ) � � 1 , π ( w i | ... ) α ( w i | x i ) = min π ( x i | ... ) – if accepted: x t +1 := x i + δ i ( α = 1 if π ( w i | .. ) > π ( x i | .. ) ) i – if rejected: x t +1 := x i i

Hierarchical Mixed Effects Model • First level ǫ ij ∼ N ( 0 , Σ 1 ) y ij = f ( ϑ i , d i , D i , t ij ) + ǫ ij , ϑ i = (log k i , logit ϕ i ) , d i = (log p bi , log x 1 b , log x 2 b ) • Second level t − Student( ν, ¯ ∼ ϑ, Σ 2 ) ⇒ ϑ i N (¯ λ i ∼ Ga( ν 2 , ν ∼ ϑ, Σ 2 /λ i ) p ( λ i ) , 2 ) N ( ¯ d i ∼ d, Σ 3 ) • Third level � ρ l , ( ρ l R l ) − 1 � Σ − 1 ∼ Wishart l = 1 , . . . , 3 l ∼ N ( c, C ) µ η ∼ N ( M, S )

Hierarchy in FFA and KB Model ρ 2 ρ R ν c C M S R 2 3 3 η Σ λ µ Σ 2 3 ρ 1 R 1 ϑ Σ X d 1 t y D

Model Validation (Example)