Differential equation methods for simulation of GFP kinetics in non-steady state experiments.

Genetically encoded fluorescent proteins, combined with fluorescence microscopy, are widely used in cell biology to collect kinetic data on intracellular trafficking. Methods for extraction of quantitative information from these data are based on the mathematics of diffusion and tracer kinetics. Current methods, although useful and powerful, depend on the assumption that the cellular system being studied is in a steady state, that is, the assumption that all the molecular concentrations and fluxes are constant for the duration of the experiment. Here, we derive new tracer kinetic analytical methods for non-steady state biological systems by constructing mechanistic nonlinear differential equation models of the underlying cell biological processes and linking them to a separate set of differential equations governing the kinetics of the fluorescent tracer. Linking the two sets of equations is based on a new application of the fundamental tracer principle of indistinguishability and, unlike current methods, supports correct dependence of tracer kinetics on cellular dynamics. This approach thus provides a general mathematical framework for applications of GFP fluorescence microscopy (including photobleaching [FRAP, FLIP] and photoactivation to frequently encountered experimental protocols involving physiological or pharmacological perturbations (e.g., growth factors, neurotransmitters, acute knockouts, inhibitors, hormones, cytokines, and metabolites) that initiate mechanistically informative intracellular transients. When a new steady state is achieved, these methods automatically reduce to classical steady state tracer kinetic analysis.

INTRODUCTION

Fluorescent probes, with their characteristic rapid response times and enormous potential dynamic range, have long been flexible and productive sensors of pivotal cell physiological and biophysical variables (Tsien, 1980; Paradiso ; Farkas ). In cell biology, the power of fluorescence imaging has emphasized the use of fluorescent tags or tracers to follow the cellular itineraries of specific molecules. Moreover, the possibility—unique to fluorescence—of switching off the tracer on command, by photobleaching has been the source of many cell biological insights. Extraction of useful information about diffusion by quantitative analysis of photobleaching recovery kinetics dates from classic work published in the late 1970s (Axelrod ; Edidin ). With the advent of green fluorescent protein (GFP) genetic constructs (Chalfie ; Heim ) it became possible to extend the determination of diffusion coefficients to intracellular membranes and compartments (Cole ; Seksek ; Partikian ). Consequently, analysis of photobleaching kinetics for studies of protein diffusion and binding has been extensively developed. Excellent methods papers and reviews are available (Axelrod ; Sprague and McNally, 2005; Beaudouin ; McNally, 2008; Mueller ; Stasevich ; Blumenthal ). Typically, experiments of this kind involve data collection on a time scale of seconds, and the principal challenge is to separate the kinetics of binding from the kinetics of diffusion.

A transition soon followed from studies of diffusion/binding to studies of endomembrane and protein trafficking on time and space scales much greater than those for which diffusion is a concern (Presley ). Once feasibility was established, a natural next step was quantitative analysis of the trafficking kinetic data (Hirschberg ). Thus, the transition was made from the mathematics of diffusion and binding to the mathematics of tracer kinetics. Closed-form solutions, so effective for simple diffusion and binding models, are impractical or even impossible as complexity increases, and therefore gave way to the ordinary differential equations of compartmental analysis.

RESULTS AND DISCUSSION

Steady state kinetic analysis

Quantitative analysis of GFP-tagged protein kinetics has generally assumed, either implicitly or explicitly, that the cells, from which the data are collected, are in a steady state (Hirschberg ). A simple definition of steady state is that nothing is changing with time. A more technical definition is that dS/dt = 0 for j = 1,…,n, where the S are the individual state variables in the biological system of interest. In a cellular system, state variables (or simply states) are the concentrations (nM) or abundances (molecules/cell) for each molecule in each cellular compartment (e.g., cytosol, endoplasmic reticulum [ER], nucleoplasm, mitochondrial matrix, plasma membrane). In any given experiment, investigators are typically interested in a subset of the full collection of states, but the size of that subset seems to increase every year as cell biology aims at more and more comprehensive models of cellular function. The assumption of steady state is often entirely reasonable. If the culture medium contains the essential substrates, degradation of all molecules is balanced by synthesis, entropy is being produced at the minimum rate, and the concentrations, fluxes, and membrane potentials are constant. If a steady state is attained, and the other requirements for classical tracer kinetics are met, then the dynamics of the tracer are guaranteed to be described by a model consisting of linear ordinary differential equations with constant coefficients. We will call these the steady state tracer differential equations:

where represents the abundance of the tracer (e.g., GFP-fusion protein) in the jth compartment (e.g., ER lumen) of the model and k is the rate constant characterizing the process that transfers the tracer molecule into compartment i from compartment j. These rate constants quantify the fraction of the molecules in compartment j that move to compartment i per unit time and have units of inverse time (e.g., min−1). Informally, a rate constant of 0.025 min−1 can be reported as 0.025 pools/min, or 2.5% per minute. These forms are entirely equivalent.

The steady state tracer differential equations 1 express conservation of mass for each tagged molecule in each cellular compartment; the first term on the right-hand side is the sum of all the fluxes (molecules min−1 cell−1) entering compartment j and the second term is the sum of all the fluxes leaving. Moreover, the underlying biochemical–cell biological system, because it is (or is assumed to be) in steady state, requires only linear algebra for its analysis (Berman ; Jacquez, 1996).

Extension to non–steady state experiments

Although steady state analysis is extremely useful, cell biological and biochemical experiments often involve perturbations of the steady state. These include changes in temperature or pH; addition of drugs, hormones, or growth factors to the medium; and changes in nutrient or ionic composition by replacement of the culture medium. Perturbations are chosen to activate or inhibit processes of interest to the investigators. They induce changes that drive the cells out of their initial steady state, through a transient whose duration is generally not known a priori, and possibly into a new steady state, if the perturbation is sustained or the system is nonlinear and bi- or multistable.

Because the underlying cell biological system is no longer in steady state, we naturally end up with two closely related systems of differential equations. And those close relations are the key to data analysis. The first system we will call the chemical system differential equations, and the second system we refer to as the non–steady state tracer differential equations.

Chemical system differential equations.

The chemical system equations are so-named to emphasize that they represent the chemical kinetics of each process in the model so that the dynamics of each state can be calculated based on changes in its inputs and outputs. It is helpful to include both native molecules and traced molecules in each chemical system state in order to account for those, not infrequent, cases where the addition of tracer molecules is substantial enough to perturb the steady state. This happens, for example, when transfection of the GFP-encoding genetic construct results in overexpression of the encoded protein, and also in modern stable isotope tracer experiments where a large dose of tracer-labeled molecules is often necessary to ensure the tracer is detectable in its target pathways.

In general, the kinetic laws or rate laws for these processes (P) are nonlinear functions, P(S(t)), where S(t) is the state vector of the system. In other words, the rate laws tell us that the fluxes (molecules/min/cell) depend on the current abundances (molecules/cell) or concentrations (nM) of at least one, and often many, of the molecules present in the cell. A simple example of such a rate law is the Michaelis-Menten equation with a competitive inhibitor (Sauro, 2012). The well-known characteristic constants, kcat, Km, and Ki quantify the relationships among the molecules that cause changes in the flux through an enzyme. Applying conservation of mass to each molecule in the modeled system, we obtain a system of nonlinear ordinary differential equations called the chemical system differential equations:

where S (in, say, molecules/cell) is the jth element of the chemical system state vector, S(t), P is the flux (in, say, molecules min−1 cell−1) into the jth state from the ith state, and n is the number of states. There is, of course, no need for the state vector, S(t), to include every molecule in every pathway your cell expresses. On the other hand, S(t) must include all the molecules whose abundances or concentrations you wish to plot, as well as all those you plan to measure or whose concentrations you plan to manipulate experimentally.

Applying the fundamental assumption of tracer kinetics, namely, indistinguishability between parent molecule (e.g., protein of interest) and tracer molecule (e.g., GFP-tagged protein of interest), we then write all the tracer fluxes, (in, say, tracer molecules min−1 cell−1), as

Notice, particularly, that the coefficient of in Eq. 3 is not a constant. Because the experiment disturbed the steady state, all n of the chemical system states, including S, are liable to be changing in time. If we had remained in a steady state, both P and S would be constant and the non–steady state tracer differential equations would revert to their classical linear constant coefficient form with as introduced above in Eq. 1. In other words, the equations for non–steady state tracer kinetics automatically reduce to those for classical tracer kinetics whenever the underlying biological system returns to a steady state.

Importantly, Eq. 3 is the critical link between the chemical system differential equations and the non–steady state tracer differential equations; it calculates each tracer flux based on the non–steady state chemical system flux, and the probability, , that a molecule in the source state is a tracer molecule. Notice that drugs and other experimental perturbations do not appear explicitly in the non–steady state tracer differential equations except insofar as their effects are conveyed by causing perturbations of chemical system states and fluxes in Eq. 3. Tracer fluxes thus change automatically and correctly in response to perturbations of the underlying chemical system. Eq. 3 is the key result of this work.

Non–steady state tracer equations.

Now applying conservation of mass to all the tracers, the non–steady state tracer differential equations become where the are the tracer fluxes from Eq. 3.

When only one tracer molecule is used, Eqs. 2, 3, and 4 provide a complete mathematical model for the dynamics of both the tracer and the underlying non–steady state cell biological system. More tracers in the same experiment are easily accommodated by including a set of Eqs. 3 and 4 for each new tracer. When multiple tracers are used in the same experiment, the modeler’s work is easier if the software provides tools to specify which states and which processes are traced by each tracer. Some of these assignments can be made computationally, but especially for processes that yield two or more products, no algorithm can make these decisions as efficiently as a human investigator.

Applying experimental protocols to the mathematical model

Equations 2, 3, and 4 are extremely general; they represent not only your current understanding (or theory) of how the biological system works, but also provide the foundation upon which any number of specific experiments can be built. Because the goal of kinetic modeling (Beard and Kushmerick, 2009; Phair, 2014) is to test a theory against the results of particular experiments by comparing model solutions to experimental data, we need an explicit procedure for modifying the model equations to reflect exactly what was done, experimentally, to the cell biological system we are studying. The simplest and least error-prone approach is to replicate, in the model equations, exactly what was done in the laboratory. In effect, this translates the time line of the experimental procedures to a time line of quantitative changes imposed on the model equations. How to do this for each of the most widely applied fluorescent tag experimental designs is the topic of the next few sections.

Adding the tracer.

When the tracer molecule is encoded by a transiently transfected genetic construct, the outcome of an experiment can depend on the strength of the promoter and the time allowed for transcription and translation. For fusion constructs driven by strong promoters, like CMV or SV40, substantial overexpression is not uncommon. Even when this level of expression does not alter your protein’s localization or function, the experiment is immediately a non–steady state perturbation. Every cellular process for which this protein was limiting will undergo an increase in its effective Vmax until a new steady state is reached or some other molecule becomes limiting.

Suppose the plasmid-encoded protein is cytosolic and its symbol in the model state vector is S. To include transcription and translation from the transfected plasmid in the model equations we modify both the chemical system differential equations and the non–steady state tracer differential equations by adding the same protein translation process to the right-hand sides:

Units for this added term are flux units (e.g., molecules min−1 cell−1). It is important to add to the chemical system equations because this accounts for possible overexpression and resultant nonlinearities such as cooperativity and saturation. Including this ongoing production term often explains unanticipated fluorescence signals late in an experiment. When ongoing production obscures a key feature of your data, photoactivation (see below) can be a useful workaround.

Modern modeling software tools support experimental perturbations at specific times and over specific intervals of time relative to a user-defined zero of time. In this case, the time of transfection would be specified. This is the significance of in Eqs. 5 and 6.

Generally, little attention is paid to the tracer distribution phase, the period between ttransfection and the start of the experiment. The usual rationale is that the fluorescence data will always be normalized to whatever value is attained at the start of the “real” experimental protocol. It is, however, good practice to record the clock times of transfection and experiment start and to keep the difference between them as constant as possible for different experimental days. The data analyst will want to know the duration of the distribution phase because must be “turned on” at some fixed time, ttransfection, relative to the start of the first perturbation (e.g., drug addition, fluorescence recovery after photobleaching [FRAP]) of the experimental protocol and changing that time can change the contribution of slow pools to the model solution. Imagine, for example, that your protein participates in the formation of two different complexes. Imagine further that one of those complexes is relatively labile and the other is so stable that FRAP studies generally report an immobile fraction. The three FRAP curves in Figure 1 are simulations of three different tracer distribution times (8 , 12 , and 24 h) for exactly the same binding and bleaching para­meters. Short tracer distribution times clearly underestimate the magnitude of the slow (“immobile”) binding pool, and variable distribution times will add unnecessary noise to your results. It is also true that different amounts of incorporated plasmid will lead to different levels of expressed protein, which can be important if changing expression changes the identity of the rate-limiting step due to nonlinearity, but in the absence of such a nonlinearity, consistency of tracer distribution time is the more important aspect of experimental design when “immobile” fractions are nonzero.

FIGURE 1:

Sensitivity of “immobile fraction” to tracer distribution time. The slow (“immobile”) binding site has a residence time of ∼22 h. The faster binding site has a residence time of 2.7 min. Tracer distribution times for the three FRAP experiments are as indicated in the figure. FRAP is initiated at t = 0. The different apparent immobile fractions of 36%, 42%, and 51% for the different tracer distribution times, despite identical binding kinetics, emphasize the importance of consistency in the time between transfection of the GFP construct and the FRAP.

Changing extracellular stimuli.

In cell biology, most non−steady state experiments are initiated by pipetting a small volume of concentrated solution into the existing culture medium or by completely removing the medium and then replacing it with a medium of different composition, sometimes with an intervening wash. Increasingly, automated microfluidics-based tools are being used to control the composition of the extracellular medium as desired (Salieb-Beugelaar ; Ly ; Haque ). No matter whether changes are effected manually or by servomotors, replicating experimental perturbations in the model is a universal requirement for analysis of non−steady state experiments.

First, any variable that is to be manipulated by the experimentalist must appear in the chemical system differential equations. The numerical value of a state representing a drug, or a cytokine, or a hormone may be, and often is, zero during the prestimulus control period, but whenever dictated by the experimental protocol the simulation software must be capable of a stepwise change to a new value, or an incremental bolus addition if the experimental stimulus was pipetted into the dish or microinjected into the cell under observation.

If the cells reside in a flow-through chamber, a sensible approach is to model the chamber and the flow especially if the perfusate is to be collected and sampled for measurements of molecules or tracers released by the cells. For any measurements made in the extracellular medium or perfusate, it is highly desirable to have an accurate count of the number of cells in the chamber so that these measurements can be integrated with those that are made on single cells.

Every model has a boundary. Processes and events outside this boundary do not appear in the model equations. The boundary thus defines the scope of the model. Quite often, however, there are several states that reside on the boundary and a few processes that cross the boundary. States on the boundary do not have their own defining differential equations; if they did, they would be in the model, not on the model boundary. Such states must be defined for the model to be solved, and this is often done by measuring the concentration represented by that state and using interpolation between the measured points to define a forcing function for the time course of the state in question. Extracellular stimuli not under the direct control of the investigator are much more common in whole animal experiments, but are a frequent feature of coculture experiments as well. A paracrine cytokine released by an unmodeled cell type may be measured and used to define a forcing function on the model boundary for the cell type being studied.

Photobleaching: FRAP.

One of the most useful features of GFP tracers is that they may readily be photobleached, effectively switching off their ability to fluoresce.

To see how the chemical system and tracer differential equations are modified to replicate a FRAP protocol, we need only consider the model state where the fluorescent tracer is synthesized and the state where it is bleached. Consider, for example, a secretory protein, synthesized in the ER and bleached in the Golgi apparatus. If S is the state representing the protein in the ER and S represents the same protein in the Golgi, then the chemical system and tracer differential equations that must be modified to reproduce the FRAP protocol are 7, 8, and 10. Equations 7 and 8 include the terms for biosynthesis of the fluorescent tracer in the ER as described in the section on Adding the tracer, above. Equation 10 includes the term for photobleaching.

Because the kinetics of transition among the various energy states in a fluorophore’s photocycle are fast relative to photobleaching kinetics, it is standard practice to treat the photocycle as a rapid equilibrium subsystem (Wustner ) and describe Pbleach as a first-order removal process characterized by a single intrinsic bleach rate constant. In this case,

By fitting both the prebleach and postbleach data to , it is a simple matter to extract from the experimental data. Once again, the ability of a software tool to start and stop a process at specific times is important to accurate replication of the experimental protocol. In this particular case, tstop is only slightly greater than tstart because the duration of the bleach pulse is typically short on the time scale of the FRAP experiment, so it is essential that the numerical integrator ensures it does not take too large a time step and miss the photobleaching event entirely.

The photobleaching term, , is added to the right-hand side (RHS) of all tracer differential equations representing traced states within the photobleached region of interest (ROI). No corresponding term is added to the RHS of the chemical system ordinary differential equations (ODEs; Eq. 9) because photobleaching does not alter the ability of the tagged molecules to participate in chemical reactions. Solving the full system of differential equations is superior to fitting only the recovery curve (the data for ) to an arbitrary sum of exponentials because it forces the recovery kinetics to be explained by the same processes and parameters that account for all the other data constraints, such as measured concentration or abundance data. Moreover, the resulting parameter values provide quantification of those key physiological processes. A sum of exponentials is a perfectly good mathematical model but its coefficients and its exponents cannot readily be assigned cell biological meaning.

Photobleaching: FLIP.

To replicate a fluorescence loss in photobleaching (FLIP) protocol, the procedure is essentially identical to that for FRAP. The differential equations to be modified are identical to Eqs.–11. The duration of the photobleach is the biggest difference; tstop will be much greater than tstart and may, in fact, extend the bleaching period to the end of the experiment. Typically, the collected and fitted data will correspond to a state or states different from the bleached state(s). When the measured ROI includes more than one tracer state, an ancillary equation such as is necessary. Here, F is the fluorescence remaining in the ROI, knorm is a normalization constant, and the sum includes all tracer states within the defined ROI.

Photoactivation.

Photoactivation is widely recommended for studies of protein turnover because there is no constant biosynthesis of the fluorescent tracer to confound the fluorescence data. Several approaches are possible for replicating a photoactivation process in the context of a differential equation model. Here, we present the one that most closely represents what actually happens experimentally and biophysically. This straightforward approach requires three sets of differential equations: 1) the chemical system differential equations with the addition of a term representing biosynthesis of the nonfluorescent photoactivatable tracer (Eqs. 12 and 15), 2) the non−steady state tracer differential equations for the invisible (not yet photoactivated) fluorescent tracer (Eqs. 13 and 16), and 3) the non−steady state tracer differential equations for the visible (photoactivated) fluorescent tracer (Eqs. 14 and 17).

Here, represents the biosynthesis of the nonfluorescent (lowercase pa) photoactivatable tracer. represents the photoactivation of the invisible tracer in the Golgi. The term is negative in Eq. 16 because photoactivation removes some flux of tracer from the invisible pool; the term is positive in Eq. 17 because photoactivation adds the same flux of tracer to the visible pool. The duration of the photoactivating laser pulse is tstop – tstart.

There are three classes of photoactivatable fluorescent proteins: photoactivated, photoconverted, and photoswitched. The same system of differential equations can be used for each of these. For the photoconversion case, the pa differential equations represent the protein fluorescing at one of its wavelengths while the PA differential equations represent the protein after it has been photoconverted to fluoresce at a different emission wavelength.

Some fluorophores will be partially photobleached in the process of photoconversion, so it may be necessary to add a photobleaching term as well as the photoconversion term. Photoactivation, photoconversion, and photoswitching all take place on a time scale of 40 s or less. If the model operates on this time scale it is prudent to consider more detailed kinetics of these processes. Otherwise, as for photobleaching, the photoactivation flux, , where is approximately linear in laser power.

As a specific example that emphasizes how photoactivation is modeled using the differential equation approach, consider the simple model of albumin synthesis and secretion from hepatocytes in Figure 2. A transfected photoactivatable albumin tracer is translated on bound ribosomes and translocated into the ER lumen. A steady state is established and then photoactivation of the Golgi ROI begins at t = 0. Figure 3 displays the simulation results for the total chemical system albumin (Figure 3A), the invisible albumin-paGFP (Figure 3B), and the visible albumin-PAGFP (Figure 3C). The appearance of albumin-PAGFP in the ER is due to retrograde traffic and will have an impact on the observed decline in the Golgi signal, especially if transport from ER to Golgi is slower than anterograde export from the Golgi. A comparison of Figure 3, B and C, emphasizes the explicit conversion of invisible to visible albumin-paGFP, which assures that the abundance of albumin-paGFP available for photoactivation in any compartment is correct even if additional photoactivations are initiated between t = 0 and t = 120 min when the first photoactivation transient will still strongly determine the response to a second photoactivation.

FIGURE 2:

Model for steady state photoactivation example. Simple model of albumin secretion from hepatocytes. Arrows represent processes of albumin translation, vesicular traffic from endoplasmic reticulum (ER) to Golgi (G), retrograde traffic from G to ER, and vesicular traffic from G to the plasma membrane for exocytosis to the extracellular fluid (ECF). The proposed computational methods are applied to this system to show how model solutions (Figure 3) are obtained when the albumin-paGFP construct is transfected at t = −1440 min, and photoactivation of Golgi albumin-paGFP is carried out at t = 0 min.

FIGURE 3:

Simulation of steady state photoactivation example. Transfection with albumin-paGFP at t = −1440 min. Photoactivation of Golgi albumin-paGFP 9 s starting at t = 0. (A) Solution of chemical system differential equations for total (endogenous + tracer) albumin in endoplasmic reticulum (ER) and Golgi. (B) Solution of tracer differential equations for invisible albumin-paGFP in ER and Golgi. (C) Solution of tracer differential equations for visible (photoactivated) albumin-PAGFP in ER and Golgi.

Tracer experiment in a non–steady state: a full example.

To cement the ideas presented in this article, it may be useful to consider a specific non−steady state example in detail. Figure 4 represents a mechanistic model of a simple generic growth factor signaling cascade and its action on a transcription factor. This model aims to be complex enough to illustrate the details of the differential equation method while at the same time including only enough processes to provide a realistic example. More realistic models of cell signaling can easily be accommodated.

FIGURE 4:

Hypothetical simplified model of non−steady state transcription factor (TF-GFP) system. Labeled vertical “swim lanes” define physiological compartments or places. Rectangles represent molecules and molecular complexes present in those places (culture medium, plasma membrane, cytosol, and nucleus). Black arrows represent the three main types of biological processes: binding, chemical reactions, and transport. Dashed green lines terminated with + adjacent to a process indicate enzymatic catalysis. A dotted green line indicates an activator. Processes with no endpoint molecule (e.g., TF degradation) and processes with no starting point molecule (e.g., TF synthesis) are crossing the model boundary. They represent degradation to unmodeled products or synthesis from unmodeled substrates, respectively.

Beginning at the left side of the model diagram, a growth factor (GF) in the culture medium binds to a plasma membrane growth factor receptor (GFR). The off-rate constant is consistent with the usual nM KD. The ligand-bound receptor (GF:GFR) catalyzes the phosphorylation of inactive protein kinase (MAPK). A phosphatase can reverse this process. Activated kinase (MAPK-Pi) phosphorylates inactive transcription factor (TF) and the resulting activated TF (TF-Pi) is transported to the nucleoplasm where it can bind to and be released from TF response elements (TFREs) in the upstream control regions of TF-responsive genes. One such gene is modeled: the TF-Pi:TFRE phosphatase, which provides self-regulating feedback by dephosphorylating and releasing the bound TF, which is transported through nuclear pores and returned to the cytosolic TF pool. Molecules and complexes in green rectangles are traced with TF-paGFP.

Because the goal of this article is to illustrate how this system is modeled using the differential equation approach, the complete set of non−steady state chemical system differential equations defining this model, plus the non−steady state tracer differential equations for the TF-paGFP photoactivation experiment are included in the Supplemental Material.

Figure 5 displays simulation results for this model. After relaxing to a steady state, a bolus of a GF is added to the culture medium at t = 0 initiating a non−steady state transient. GF immediately increases, then decreases as it is bound to plasma membrane GFR and internalized; free GFR decrease accordingly (Figure 5A). Bound GFR (Figure 5A) increase and drive conversion of cytosolic inactive MAPK (Figure 5B) to its phosphorylated active form (Figure 5B). This active kinase then phosphorylates a cytosolic TF and TF-Pi (Figure 5C) increases in the nucleus where it also forms a complex (TF-Pi:TFRE; Figure 5C) with TF TFREs. Among the many possible genes activated by this complex, the model emphasizes a phosphatase that dephosphorylates TF-Pi to TF (Figure 5C) and terminates TFRE binding, providing self-limiting negative feedback. Free TF is returned to the cytosolic TF pool. In Figure 5D this GF experiment is run three times in a system expressing photoactivatable GFP-labeled TF (TF-paGFP). Each experiment includes photoactivation of nuclear TF-paGFP at a different time. Notice that the kinetics of photoactivated TF-PAGFP are qualitatively and quantitatively different when photoactivation is initiated at t = −100 in the pre-GF steady state compared with the same photoactivation at t = 20 during the early GF-induced transient, and different yet again when the same photoactivation is initiated later in the GF-initiated transient at t = 100 min (Figure 5D). These kinetic differences would not be predicted by a model that does not include quantitative representation of the underlying cellular feedback systems. Experimental data on these kinetic differences thus contain quantitative mechanistic information about regulatory nuclear processes that can be extracted by fitting the model solutions to the data (Phair, 2014).

FIGURE 5:

Simulation of nonlinear non−steady state signaling cascade/transcription factor model shown in Figure 4. The simulated experiment is a response to growth factor exposure at t = 0 plus photoactivation of PAGFP-labeled transcription factor before and at two times after growth factor addition (see details in the text). Full model equations are shown in the Supplemental Material. GF, growth factor; GFR, GF receptor; PM, plasma membrane; MAPK, mitogen activated protein kinase; MAPK-Pi, phosphorylated (activated) MAPK; TF, transcription factor; TF-Pi, phosphorylated (activated) TF; TFRE, TF response element.

Data analysis

Although not the focus of this article, it may be helpful to list briefly the software tools one would need to carry out this analysis for real experiments. First, a solver or numerical integrator is required. There are free tools supported by academic laboratories such as Virtual Cell, Berkeley Madonna, WinSAAM, and Copasi, as well as commercial products such as Matlab, Mathematica, SAAM2, and ProcessDB. Different software tools have different user interfaces and support different specialized functions. All of them include one or more numerical differential equation solvers, and all of them provide fast solutions based on state-of-the-art algorithms. Moreover, given the same equations, all will output the same answers. The reason there are so many choices is that scientists come to biological modeling with many different backgrounds and every developer has in mind a specific user group and a particular modeling philosophy. There is not, and may never be, a single universal approach to cell biological modeling.

The next question is, What, in terms of model variables, was measured experimentally? Each software tool must have a means of adding ancillary equations to the differential equation model. The most important of these are the equations that calculate observed variables from model variables. For example, in the TF model described above, measurements in the nuclear ROI will not be able to distinguish TF-PAGFP from TF-Pi-PAGFP. Indeed, the data from the nuclear ROI will have to be fitted to the sum of all three TF-PAGFP species in the nucleus. Many other relationships between model variables and observed variables are found in published models.

Ancillary equations have other valuable applications. They make predictions for variables we cannot yet measure, and they help design experiments to measure them. They do unit conversions from model variables to related variables that readers want to see. They do normalizations and they calculate figures of merit that often serve as useful constraints.

The final software tool is optimization. Here there is enormous variety in what is offered by different tools. Experience suggests choosing a tool with both a classical local gradient optimizer, and one of the many global optimizers. The goal of parameter optimization is to find a set of numerical values for the parameters of a model that succeeds in minimizing an objective function that depends on the sum of the squares of differences between model solution and experimental data or maximizing the likelihood of the experimental data given the model and parameter values. Successful gradient optimization provides statistical information on mechanistic parameters of biological interest, and global optimization has the potential to treat the model as a hypothesis and test it by asking whether there is any combination of numerical parameter values that can satisfactorily account for the experimental data.

Finally, it may be worth repeating that the differential equation procedures detailed above for modeling the addition of tracers, FRAP, FLIP, and photoactivation in non−steady state experiments are equally applicable to steady state experiments. Consequently, the modeling approach described here is general enough for the analysis of kinetic data from a wide variety of steady and non−steady state experimental protocols involving fluorescent tracers. The equations derived in this article and in the Supplemental Material should allow interested investigators to implement these techniques in their own research environments.

Click here for additional data file.