Mathematical Modeling of Cellular Signaling in Macrophages: Lipid Signaling Kinetics

Mary Ann Horn National Science Foundation & Vanderbilt University Joint work with Hannah L. Callender, Department of Mathematics & H. Alex Brown & the Brown Lab, Department of Pharmacology

Collaboration with Alliance for Cellular Signaling

• Alliance for Cellular Signaling (AfCS): “Determine how cells process signaling information in a context dependent manner” • AfCS: Screen a large number of ligands (signaling molecules) and analyze their effects on the RAW 264.7 cell (a macrophage-like cell line derived from mice; a macrophage is a type of white blood cell that surrounds and kills microorganisms, removes dead cells, and stimulates the action of other immune system cells). • Focus: Construction of a comprehensive mathematical model for the signaling pathway of one ligand, uridine 5’-diphosphate (UDP), in the RAW 264.7 cell.

– UDP is a common signaling nucleotide often released when cells undergo some type of trauma or mechanical stimulus • Modeling approach: Utilize a system of nonlinear ordinary differential equations designed to model the time-dependent behavior of key players in the pathway.

Basics of Cellular Signaling

Signaling molecules bind to cell surface receptors  Initiate a series of intracellular reactions  These reactions regulate virtually all aspects of cell behavior (including metabolism, movement, proliferation, survival, and differentiation) • Breakdown in the signaling pathways that control normal cell proliferation and survival result in many types of cancer  Need for better understanding  Mathematical modeling can provide insight into these signaling pathways

Canonical UDP signaling pathway

1. Ligand, or stimulus (UDP), binds to a specific G-protein coupled receptor (P2Y

6

).

3. G



q ·

GTP activates phospholipase C  3 (PLC-  3).

GDP GTP

2. The ligand bound receptor causes the exchange of GDP for GTP on the G 

q

subunit of the G-protein.

Ca 2+ Ca 2+ GTP Ca 2+ IP 3 IP 3 Ca 2+ Ca 2+ Ca 2+ Ca 2+ Ca 2+ 5. IP 3

binds to a specific receptor on releasing sequestered Ca

2+

.

Ca 2+ Ca 2+

4. Active PLC-  3 cleaves phosphatidylinositol 4,5 bisphosphate (PIP

2

) into (IP

3

) and diacylglycerol

Mathematical Model

• The model consists of a base set of 10 coupled nonlinear ordinary differential equations (ODEs).

• The ODEs are constructed using mass action kinetics to approximate diffusion, binding and unbinding of molecular species in various compartments of the RAW cell.

 effective rates of production and degradation are proportional to the concentrations (or number of molecules) of participating reactants • Four modules: – Receptor regulation – G-protein casade – Species-specific DAG dynamics – Cytosolic Ca 2+ dynamics

Model Equations: Receptor regulation

• Receptor-ligand binding (based on work by Lemon 2003): – Number of “activated” surface receptors at time t: – Number of “inactivated” surface receptors at time t: – Concentration of ligand (constant):

Nondimensionalized equations

Model Equations: G-protein activation

• G-protein activation: – Number of activated G-proteins (G  GTP) at time t:

Nondimensionalized equations

Model Equations: PIP

2

hydrolysis and replenishment

• PIP

2

hydrolysis and replenishment:

GTP GTP Ca 2+ IP 3

– Number of PIP 2 molecules at time t: – Initial number of “free” PIP 2 :

Nondimensionalized equations

Model Equations: IP

3

production and degradation

• IP

3

: • Concentration of IP 3 at time t (in  M):

Nondimensionalized equations

Model Equations: DAG

• DAG (separate ODEs for each species considered): • Concentration of DAG at time t:

Nondimensionalized equations

Model Equations: Li and Rinzel (1994) Ca

2+

• Free cytosolic Ca

2+ :

• h = fraction of IP 3 channels not yet inactivated by Ca 2+ : • c

0

= concentration of total free Ca 2+ per cytosolic volume:

Timeseries data for the model

• Data from collaborators at UTSW (Sternweis Lab) – Inositol 1,4,5-trisphosphate (IP

3

) production post 25  M UDP stimulation in the RAW cell – Cytosolic calcium (Ca

2+

) release post stimulation with 25  M UDP in the RAW cell.

• Diacylglycerol (DAG) data H. Callender collected from a novel method of quantitative analysis of multiple species of DAG, developed in the Brown Lab. Reference: Callender, H. L., et al. Quantification of Diacylglycerol Species from Cellular Extracts by Electrospray Ionization Mass Spectrometry Using a Linear Regression Algorithm. Anal. Chem. 79 (2007) 263-272.

A.

IP

3

and Ca

2+

data from AfCS

• IP  3 25 response (in pmoles per 100  M UDP • Points represent the average of four experiments.

• Cytosolic calcium response (in  M) in the RAW 264.7 cell to 25  M UDP. • The graph displays 43 experimental repeats.

B.

What is DAG?

• DAG is a cellular second messenger molecule which plays an important role in initiating various changes in cell behavior, including cell activation, differentiation, proliferation and tumor promotion. • There are many different species of DAG, depending on the number of carbons and number of double bonds in the fatty acyl chains, and different species can have different cellular functions.

• Example:

32:0 DAG

O O O O Two fatty acyl (hydrocarbon) chains in every DAG species Each corner represents carbon OH

32:0 DAG: 18 carbons in first chain, 14 in the second, no double bonds in either chain

Kinetics of mono/di unsaturated DAGs

• Time based behavior of four mono/di unsaturated DAG species after addition of 25  M (solid red squares) and 0.25  M (solid green triangles) UDP • Time points contain a minimum of nine replicates performed on three different experimental days

32:1 DAG 34:1 DAG - 25



M UDP - 0.25



M UDP - 25



M UDP - 0.25



M UDP 0 5 10 15 20 Time (min.) 25 30 - 25



M UDP - 0.25



M UDP 34:1 DAGe/p 0 5 10 15 20 Time (min.) 25 30 - 25



M UDP - 0.25



M UDP 36:2 DAG 0 5 10 15 Time (min.) 20 25 30 0 5 10 15 20 Time (min.) 25 30

• • •

Kinetics of PUFA containing DAGs

Time based behavior of four polyunsaturated fatty acid (PUFA) containing DAG species after addition of 25  M (solid red squares) and 0.25  M (solid green triangles) UDP Time points contain a minimum of nine replicates performed on three different experimental days Mono/di unsaturated DAG species give a larger increase than polyunsaturated fatty acid (PUFA) containing DAGs

38:5 DAG 38:3 DAG - 25



M UDP - 0.25



M UDP - 25



M UDP - 0.25



M UDP 0 5 10 15 Time (min.) 20 - 25



M UDP - 0.25



M UDP 38:4 DAG 25 30 0 5 - 25



M UDP - 0.25



M UDP 10 15 Time (min.) 20 36:4 DAG 25 30 0 5 10 15 Time (min.) 20 25 30 0 5 10 15 20 Time (min.) 25 30

Differential DAG Kinetics

• Time based behavior of three DAG species with varying degrees of unsaturation after addition of

25



M

(solid red triangles) and

0.25



M

(solid green squares) UDP • Time points contain nine replicates performed on three different experimental days

34:1 DAG

Best fit with current model structure

Overall Objectives of Modeling Effort

• Predict quantitative changes in lipid species after stimulation by various ligands and ligand concentrations in the RAW 264.7 macrophage • Comparison and refinement of model output with AfCS IP 3 measurements and Ca 2+ traces as well as DAG data generated in the Brown lab • Predict in silico effects such as the effect of knock-downs, etc., on given pathway • Suggest modifications to current pathway structures

Modifications to the model

• Include an additional branch in the pathway for a second pool of DAG • Simplify Ca 2+ equations for mathematical analysis purposes

Proposed Pathway

• Note: We measure total cellular DAG levels • Initial production of DAG from the hydrolysis of PIP 2 in pool 1 (plasma membrane) is offset by phosphorylation of DAG by DAG kinase in pool 2 (Endoplasmic Reticulum? Nucleus?) to aid in the PI replacement pathway

PIP 2 POOL 1 DGK LPP PE POOL 2 IP 3 IP 4 IP 2 LPP

• Second wave of DAG is a result of resynthesis of PIP 2 , which is then hydrolyzed to form DAG and IP 3 .

Model Equations: DAG (pool 1)

• DAG (pool 1) (separate ODEs for each species considered): • Concentration of DAG from pool 1 at time t:

Nondimensionalized equations

Model Equations: DAG (pool 2)

• DAG (pool 2) (separate ODEs for each species considered): • Concentration of pool 2 DAG molecules at time t: • Baseline concentration of pool 2 DAG:

Nondimensionalized equations

Simplified Ca

2+

Equations

• Ca

2+

module (to match experimental AfCS

trace) where Ca 2+ response to 25uM UDP 0.14

0.12

0.1

0.08

0.06

0.04

0.02

0 0 100 200 300 time (s) 400 500 600

Theoretical Analysis

• Existence and Uniqueness • Positivity and Boundedness (for biological relevance) • Analysis of steady state behavior

Existence of Solutions (full model)

First we write our system of ODEs in the form:

(1) (2) (1)

Uniqueness of Solutions

Next, we denote a solution of

(1)

by , with initial condition:

(3) (1)

Since our system satisfies the hypotheses of Theorems 1 and 2 on our set of interest (for all positive time and on a positive bounded set in space), we know there exists a unique (local) solution (i.e., on some finite time interval, possibly small).

Question: Do the solutions remain positive and bounded (for biological significance)?

(3)

.

Positivity and boundedness of solutions

• We first use Theorem 3 to show positivity and boundedness of x 1 and x 2 .

• Next we use these results and the Fundamental Theorem of Calculus to show positivity and boundedness for the remaining equations. • This then ensures a global solution:

Numerical Analysis

• Parameter Estimation in SIMULINK • Sensitivity Analysis

Parameter Estimation in SIMULINK (for FULL model)

• Total number of model parameters (for full model with Li and Rinzel Ca 2+ module) = 34 – From the literature = 20 – Estimated = 14 • Receptor module: Total = 7 – From literature = 6 – Estimated = 1 (k p : rate of receptor phosphorylation) • G-protein cascade: Total = 9 – From literature = 6 – Estimated = 3 (k hyd , k rep , k d3 ) • DAG kinetics: Total = 5 (for each DAG species considered) – From literature = 0 – Estimated = 5 (all DAG parameters) • Ca 2+ module: Total = 13 – From literature = 8 – Estimated = 5

SIMULINK Details

• Unknown rate parameters were estimated using SIMULINK – Minimizes a user-specified cost function via a user-specified optimization method – Nonlinear least squares optimization method of Levenberg Marquardt was used to minimize a sum of squared errors cost function of the empirical observations and model predictions for IP 3 , Ca 2+ , and multiple species of DAG.

• Note: Although the Gauss-Newton method is often more efficient, the method of Levenberg-Marquardt has proved to be more robust.

Measured vs Simulated: 38:4 DAG and 34:1 DAG

• Timecourse of stimulation with 25  M UDP in RAW 264.7 cells.

• Solid black lines represent model simulations. • (

a

) 38:4 DAG response (representative of the response of most poly unsaturated fatty acid-containing DAG species). • (

b

) 34:1 DAG response (representative of the response of most mono- and di-unsaturated fatty acid-containing DAG species). • Data points contain nine replicates performed on three different experimental days, with error bars = 1 SEM. Units are total change in ng over baseline levels in ~8x10 6 cells.

Measured vs Simulated: IP

3

and Ca

2+

• Timecourse of stimulation with 25  M UDP in RAW 264.7 cells.

• Solid black lines represent model simulations.

• (

c

) IP 3 response in pmols per ~3.5x10

5 cells. • Points in (c) represent the average of four experiments, and error bars are 1 SEM. • (

d

) Ca 2+ response in ligand screen.

 M. Red curve is a representative Ca 2+ trace taken from the UDP experiments within the AfCS single

Model Simulations: P2Y

6

, G



GTP, PIP

2

P2Y 6 from 25  M UDP G  GTP from 25  M UDP PIP 2 from 25  M UDP • (

a

) = total P2Y 6 activated (solid line) and inactivated (dashed line) surface receptors • (

b

) = total G  GTP • (

c

) = total PIP 2 hydrolysis available for

Sensitivity Analysis

• Sensitivity analysis techniques are valuable tools designed to answer questions regarding which of the uncertain input variables is more important in determining the uncertainty in our output. • Likewise, sensitivity analysis can provide insight into which parameter should be studied in more detail in order to reduce the most variance in the model output. • The ability to answer these types of questions could lead to important insight into the design of new experiments and in determining which experiments would give us the most valuable information.

Sampling Method

• Generate a random sample of our space of input variables over a ten percent variation from each parameter's nominal value using the Latin Hypercube Sampling (LHS) method.

• Uses Standardized Regression Coefficients (SRCs) obtained by performing multiple linear regression analysis – offers a measure of sensitivity that is multi-dimensionally averaged over the entire space of parameter values. – SRCs give insight into degree of nonlinearity in the model

Computing R

y 2

values

• SRCs are only reliable measures of sensitivity when degree of nonlinearity is “small”.

• Use model coefficients of determination, R y 2 , given by • where y i is the estimate of y i obtained from the regression model.

• R y 2 ≥ 0.7 ensures SRCs are good sensitivity measures.

SRCs (as sensitivity measures) for DAG

SRCs (as sensitivity measures) for IP

3

A.

Effects of 3 most sensitive DAG parameters

B.

A. k dp2 : degradation of pool 2 DAG B. k dp1 : degradation of pool 1 DAG C. k ap2 : production of pool 2 DAG

C.

Conclusions

• We have developed a model of the UDP signaling pathway in RAW 264.7 macrophages which can predict the responses of multiple species of DAG as well as the responses of IP 3 , Ca 2+ , receptor dynamics, G-protein activation, and PIP • Full model analysis: 2 hydrolysis.

• Simplified model results: – We have obtained global existence, uniqueness, positivity, and boundedness of solutions – We have proven global stability of a unique steady state within our region of interest.

– Using SIMULINK we have estimated unknown rate parameters to obtain best fits to multiple DAG traces, IP 3 , and Ca 2+ , all in response to 25  M UDP.

– We have performed sensitivity analysis using the Latin Hypercube sampling technique in combination with standardized regression coefficients to determine which model parameters are responsible for most of the model output uncertainty.

Future Directions

• We have conducted multiple experiments to pharmacologically inhibit several different enzymes we believe to play a role in this signaling pathway, as suggested by current known metabolic pathways and by modeling results.

• The next step is to perform gene knockdowns on specific enzymes to verify results of inhibitor data.

• The model output also suggests a time delay from receptor activation to PIP 2 hydrolysis. An upcoming task is to investigate the outcomes of adding such a delay term.

• The model could be greatly enhanced by incorporating spatial dynamics, so this will also be a major focus for future research.

Acknowledgements

Mathematics Department

• Hannah L. Callender, Ph.D.

Brown Lab

• H. Alex Brown, Ph.D.

• Jeffrey S. Forrester, Ph.D.

• Mark Byrne, Ph.D.

• Anita Preininger, Ph.D.

• Michelle Armstrong • Andrew Goodman • Pavlina Ivanova, Ph.D.

• Steve Milne, Ph.D.

Collaborators

• Alliance for Cellular Signaling • UT Southwestern – Paul Sternweis, Ph.D.

– Dianne DeCamp, Ph.D.