Model-based control of the temporal patterns of intracellular signaling in silico.

The functions of intracellular signal transduction systems are determined by the temporal behavior of intracellular molecules and their interactions. Of the many dynamical properties of the system, the relationship between the dynamics of upstream molecules and downstream molecules is particularly important. A useful tool in understanding this relationship is a methodology to control the dynamics of intracellular molecules with an extracellular stimulus. However, this is a difficult task because the relationship between the levels of upstream molecules and those of downstream molecules is often not only stochastic, but also time-inhomogeneous, nonlinear, and not one-to-one. In this paper, we present an easy-to-implement model-based control method that makes the target downstream molecule to trace a desired time course by changing the concentration of a controllable upstream molecule. Our method uses predictions from Monte Carlo simulations of the model to decide the strength of the stimulus, while using a particle-based approach to make inferences regarding unobservable states. We applied our method to in silico control problems of insulin-dependent AKT pathway model and EGF-dependent Akt pathway model with system noise. We show that our method can robustly control the dynamics of the intracellular molecules against unknown system noise of various strengths, even in the absence of complete knowledge of the true model of the target system.

It is generally accepted that the functions of biological systems are determined by the temporal behavior of intracellular molecules and their interactions [1]. However, the system governing signal transduction that can be affected by extracellular stimuli is, often highly nonlinear, and it is extremely difficult to infer the molecular control system just from dose-response curves of particular intracellular molecules.

The reasons for this difficulty are multifold. These highly nonlinear systems may have multiple meta-static states and even show chaotic behaviors [2]. Dose-response relationships in such systems can also exhibit hysteresis [3]. In many cases, it is practically impossible to stabilize these systems at a desired state by simply maintaining constant the level of some controllable upstream molecule. This challenge has motivated many researchers in the field of Systems Biology to systematically examine the behavior of downstream molecules in response to various stimuli invoked by upstream molecules. These studies gave rise to the notion of temporal coding, an idea suggesting that the biological system is transferring information downstream not just via the state of the system but also via the temporal behavior of the system itself. Studies have revealed that there are patterns associated with such information transmission [1,4-8]. For biological systems, the relationship between the levels of upstream and downstream molecules is often time-inhomogeneous, nonlinear, and not one-to-one. This requires the study of the dynamics of internal states of the system that mediate the interaction between upstream and downstream dynamics. Recent studies tackled this problem and aimed to control molecular dynamics using closed-loop [9-20] and open-loop [21,22] control methods.

While both open-loop and closed-loop methods mentioned above have the common goal of causing the time course of a target molecule to follow a desired temporal pattern, they differ much in their philosophies, as they target different situations of a biological system. Open-loop methods conclude all computations in silico hoping that the response of the real biological system to the proposed upstream stimulus pattern does not diverge greatly from that of the model. Closed-loop methods on the other hand, take countermeasures against possible divergences of the real system’s response to the stimulus from predicted responses, and feedback the realtime observation of the target system’s response to the control input. In general, as suggested by studies by Milias-Argeitis et al. and Uhlendorf et al. [10,12], closed-loop methods tend to show better performance than open-loop methods.

For closed-loop control of biological systems, previous studies employed Proportional-Integral (PI) control [9,11, 14–18,20] and model predictive control [10-12,18-20]. PI control is a model-free control method, while model predictive control utilizes a mathematical model of the target system [23]. In model predictive control, previous studies used a Kalman filter to estimate the state of a linear system [10,11, 18], or a Monte Carlo filter (particle filter) to estimate the state of a nonlinear system [24]. This is helpful in attempting to control biological systems, because nonlinear mathematical models are often more appropriate to describe the dynamics underlying a variety of biological systems.

Dynamics of cellular systems are intrinsically stochastic, due to system noise often from unknown sources. Additionally, the initial concentration of the molecules of interest may vary across cells, even if they have the same genetic background [25-27]. For these reasons, control of the time course of downstream molecules tends to be much more difficult in biological systems than in industrial systems. Closed-loop methods can be effective in controlling nonlinear and biological systems because they fine-tune the control input based on real-time observable responses of downstream molecules.

In this study, we present a model-based control method that aims to make a target downstream molecule to trace a desired time course by changing the concentration of a controllable upstream molecule. The method requires a model that mimics well the target biological system. The algorithm uses predictions from Monte Carlo simulations of the model based on a particle filter algorithm, thus making it a nonlinear model predictive control method [23]. To cope with uncertainty in molecular concentrations over the cells, we initiate the algorithm with particles of different concentrations.

We applied our method to in silico control problems in an insulin-dependent AKT pathway model [5] and an Epidermal Growth Factor (EGF)-dependent Akt pathway model [4]. To mimic the real biological system, we used a stochastic differential equation with unknown system noise for the target system to be controlled. The target system differs from the deterministic model that we used in our simulation for prediction purposes. We compared our method to a standard feedback control method, i.e., a PI control, using the same set of in silico control problems.

Methods

We applied our model-based control method to several in silico control experiments. In each experiment, we attempted to make the concentration of a target molecule to trace a desired time course in the target biological system by regulating the time course of a controllable upstream molecule. We assumed that the target biological system is governed by a stochastic differential equation (SDE) with system noise of unknown origin. Based on previous studies [25-27] suggesting the presence of high variation in protein concentrations across a cell population, we also assumed that the initial concentrations of biological molecules greatly vary between cells.

Problem statement

Consider a biological system consisting of D molecular species over the time range of [0, T], and let x(t) be a D dimensional vector with x(t) representing the concentration of molecule i at time t. We would simply use x and x(·) to denote the temporal pattern x[0, T]. Here, the system is discretized with respect to time and T is an integer, so x[0, T] is a T+1 dimensional vector, and x[0, T) is a T dimensional vector. There are three groups of molecular species: (1) extracellular molecules that we can directly manipulate (stim), (2) intracellular molecules that we want to control (targ), and (3) intracellular molecules that are third party molecules affected by the first group and dynamically interacting with the second group within cell(s) (intra). We would denote them by xstim, xtarg, and xintra, respectively, and following the nomenclature we introduced above, by i=stim, targ, intra. Intra species are not the subject of our control problem. We note that intra is usually a set of multiple species, and might even be a set of virtual species that do not necessarily correspond to any actual molecular species.

Let R(xstim, xtarg,intra(0)) denote the time course of the response of x to the temporal stimulus xstim in a cell with the initial condition (xtarg(0), xintra(0)). This is a function that maps the time course xstim and the initial condition of the system to the time course of molecule i. We will omit xtarg,intra(0) and use R(xstim) when it is obvious by the context. Our control objective was to make Rtarg(xstim) close to a given reference pattern by choosing an appropriate input time course xstim. More specifically, let {t0, ..., t} be the evenly spaced time sequence at which we can both observe the concentration of xtarg and change the concentration of xstim without any delay. We denote the common inter-observation interval t+1−t by Δ. Defining the distance (i.e. error) between two time-series, y1 and y2, over the time interval [a, b] by

our goal was to find

whose search space is over every temporal pattern of input stimuli that we could realize and apply to the target biological system online. In numerical experiments, the argmin was taken by selecting the best Monte-Carlo path (particle) from the simulated set with the “mock model” that was used to predict the response of the system to the input. More details are provided in the control method section.

State-space modeling

Our model-based control method is based on state-space modeling of biological systems. In the discrete time domain, a target biological system subject to our control can be described as:

where x=(xstim, xtarg, xintra) is a state variable, f is a discretized system equation, and v is a system noise. With the assumption that there is no feedback to stim from targ or from intra, and that the system noises for targ and intra are mutually independent, we have

With this system model, our observation is given by

where y is an observation variable, h is an observation function, and w is an observation noise. As explained below, we need to include this observation function in order to infer the state of the system from possibly noisy observations. Although we do not usually know the system or the observation equation of the target system, we assume that we at least have access to their approximations. Our system model (equations (1) and (2)) and our observation model (equation (3)) represent such approximations, which, we presume, perform relatively well in predicting the behavior of the true biological system and in inferring its unobservable internal states. Specifically, we would use the system model to predict the future behavior of the target biological system and use the observation model to estimate the unobservable internal states of the target biological system based on observations. Additionally, our assumption of mutually independent system noise for intra and targ does not undermine the performance of our algorithm. This can indeed be another source of bias. However, as we show in the results, our algorithm is quite robust against the bias introduced from making wrong assumptions about the noise.

Model-based control method

Our model-based control method was based on a particle filter algorithm [28,29] (refer to Supplementary Text S1). In our in silico control experiments, we assumed the intracellular concentration of the target molecule xtarg to be observable online as , and those of other intracellular molecules to be unobservable at all times. Our model-based control algorithm consists of two steps. We iterated these two steps along the time course of the reference temporal pattern. In the prediction and control step of the simulation, we applied a short-term input xstim(t) over a time interval of length Δ to the set of possible states of the target biological system, and predicted the future behavior of the target system using the discretized approximation model we mentioned above. We will refer to the approximation model as a mock system henceforth. A possible mock system is a noiseless ordinary differential equation (ODE) system. In the numerical experiments we conducted in this study, our choice of the mock system was indeed a noiseless ODE system modeled to approximate the target SDE system. We then selected the input stimulus that the mock system predictions deemed best in making the target molecule to follow the reference temporal pattern of the target molecule. In the second filtering step, we produced plausible candidate inferences (particles) for the system state based on the predicted response of the system to the input stimulus selected in the prediction and control step (Fig. 1).

Figure 1

Schematic representation of the model-based control algorithm. (A) We search for the optimal sequence of input stimuli (green line) to make the observed values (green circles) to follow the reference temporal pattern (red line). (B) At time t4, we randomly generate (n=1, ..., N), and predict the responses (n=1, ..., N) to these stimuli based on the mock system and the current inference of xintra(t4) and xtarg(t4) (black circles). (C) In (B), we see that achieves the best response so we adopt the value for the stimulus to be applied to the target system over the interval [t4, t5). We also simulate xtarg(t5) and xintra(t5) under with the mock system (blue circles). (D) Upon the observation of , we use a particle filter to obtain the approximate empirical distributions of xtarg(t5) and xintra(t5) (black circles). We repeat the same procedure online. R̃ denotes the response of the mock system, not of the target system. N indicates the number of particles.

Here, we provide more details for each step. At the onset of the prediction and control step at time t, we generated a set of candidates for the next control input xstim(t) by adding Gaussian noise ɛ(t) of variance to the previous and actually selected control input xstim(t−1):

Equation (4) is an instance of equation (1), and the user can choose a different control equation if necessary. We adopted Gaussian random walk, ɛ for our choice of v. We assigned to each particle in the particle pool a different control stimulus from the set of candidates prepared above. We assumed that we could apply the control input continuously over each time interval, i.e., from t to t+1. It should be noted that we could change xstim only at the beginning of an inter-observation interval and that xstim was assumed to remain constant over each interval. We therefore use xstim(t) and xstim[t, t+1) alternatively in the notation. At the initial time t0, we prepared a set of control inputs xstim(t0) based on an initial distribution , which can be arbitrarily set by the user. We then simulated the aforementioned mock system

forward in time from t and for each particle with the assigned control input. Here, ϕtarg and ϕintra represent the dynamics of target molecules and intracellular molecules in the mock system, respectively. In actual implementation, the dynamical system (5) is simulated with a time interval that is much smaller than the length of the inter-observation interval, Δ. It should be noted that equation (5) behaves differently with different choices of extracellular input xstim(t). We generated an initial set of particles from a predefined initial distribution:

This initial set of particles represents the variability in the basal concentrations of molecular species across the cell population. In summary, we used our mock system model to predict the behavior of the target system at the next observation time t+1 in response to the control stimulus xstim(t).

We then selected from all the xstim(t) candidates the one that yielded the best xtarg time course on the predictive simulations based on the mock system. This approach shares much with the existing receding horizon strategy [30] and model predictive control [23]. In fact, in our full experiment, we conducted a multistep version of what we described above. Specifically, at each t, we generated short-term multistep time courses for the control stimuli xstim using the model (4) (Fig. 1B, left) up to t+−1 (τ≥1), and assigned each of them to a different particle in the particle pool. By simulating the mock system of equation (5) with each particle-xstim time course pair, we produced a set of predictions for the time course of xtarg and xintra (Fig. 1B, right) up to t+. We selected the xstim candidate whose corresponding prediction of the xtarg time course was the closest to the desired reference time course over the set of time points t+1, ..., t+. We used the squared error distance as a measure of “closeness” between a pair of time courses. We then applied the selected best control input xstim(t) to the target biological system (Fig. 1C left). In this way, in our prediction and control step, we used the measure of “closeness” to select from the randomly generated input time courses the one best realizing the time course of the target molecule that was closest to the reference temporal pattern. We then used this selected input as the optimal input.

In the next filtering step, we estimated the internal state of the target biological system based on the response of the system to the input control stimulus selected in the prediction and control step described above. After applying the selected xstim(t) to the actual system, we obtained the observation at time t+1. We would then perform another set of simulations in order to estimate the unobservable state of the biological system. With each particle (xtarg(t), xintra(t)) in the particle pool and the selected xstim(t), we simulated the mock system (5) forward in time and obtained another set of particles (xtarg(t+1), xintra(t+1)) (Fig. 1C, right). We evaluated the likelihood of the new particles, based on the observation model

(Fig. 1D). If ɛ represents Gaussian noise with variance , the likelihood is given by

We could then normalize these likelihoods and construct a probability distribution over the set of new particles. It should be noted here that the level of the observation noise ɛ does not necessarily correspond to that of the real observation noise. The variance of this Gaussian observation process, , can be arbitrarily set by the user. We then created another set of particles by resampling from the newest set of particles based on the probability distribution described above. When the initial concentration of the target molecule, , was available, as in many control situations, the filtering step preceded the prediction and control step. We used this setting later in our in silico experiments. Further details of the algorithm are described in Supplementary Text S2.

In silico application

In our in silico control experiments, we represented the real target biological system as a stochastic differential equation, that is, the noisy version of the ODE model proposed previously [4,5], since real biological systems are always perturbed by unidentified system noise. We assumed that the discretized version of the ODE, or the mock system, is known to us prior to the experiment. Thus, the mock system that we adopted in our in silico experiment coincided with the drift/ deterministic component of the target SDE system. However, in the numerical experiment, the exact initial condition of the intracellular molecules in the target system remained unknown. Because of such uncertainty associated with the target system, some discrepancy with our prediction based on the mock system was inevitable.

Application to the insulin-dependent AKT pathway in silico

Kubota et al. previously presented an ODE model that represents the dynamics of 12 molecular species in the insulin-dependent AKT pathway [5] (refer to Supplementary Text S4 for details). The model includes AKT, ribosomal protein S6 kinase (S6K), glycogen synthase kinase-3β (GSK3β), and glucose-6-phosphatase (G6Pase), which are all regulated by the extracellular insulin input (Fig. 2A). It is suggested in Kubota’s work that S6K, GSK3β, and G6Pase are selectively regulated by specific sets of patterns of AKT signal, which is in turn controlled by extracellular insulin stimuli. Here, we examined if our model-based control method could produce given reference patterns for four kinds of target molecules: AKT, S6K, GSK3β, and G6Pase, using the temporal pattern of insulin. In other words, stim in our in silico experiment was insulin, and targ was AKT, S6K, GSK3β or G6Pase. While targ represents one particular target species, xintra represents the vector of concentrations of other intracellular molecules. We represented the true target system with the SDE model

Figure 2

Control of the insulin and EGF signal transduction pathways. (A) Schematic overview of the mathematical model of the insulin-dependent AKT pathway. Red-colored molecules represent the intracellular signaling molecules to be controlled. The region surrounded by the green line constitutes the intracellular signaling network. (B) Results of the in silico control (σ=0.025) of the insulin-dependent AKT pathway. τ indicates the future value length used in the prediction and control step in the model-based control algorithm. Red trajectories show the reference temporal patterns, and blue trajectories the responses of the in silico system. Trajectories of 50 control trials are shown. The trajectories vary because of noise and the variation in the initial state of the system. (C) Schematic overview of the mathematical model of the EGF-dependent Akt pathway. (D) Results of the in silico control (σ=0.025) of the EGF-dependent Akt pathway.

where f is the ODE model presented by Kubota et al. [5], X=xintra or xtarg, and W represents independent Brownian motion with strength σ=0.0, 0.025 (default value) or 0.05. The value of σ is not very important, because in reality it is just an unknown parameter. The initial state of the in silico biological system, i.e., in silico cell, was randomly taken from a log-normal distribution with coefficient of variation of 0.25, based on a previous study [25]. It is known that protein concentration in mammalian clonal cells follows an almost log-normal distribution, with a coefficient of variation of ~0.1 to ~0.6 [26]. We also chose the parameters of the log-normal distribution so that its mean is the same as the initial value in the mathematical model presented by Kubota et al. [5]. Here, we would like to mention again that our objective was to obtain a time course of insulin, as xinsulin, that would produce the desired time course of target molecular species in the target system (8). We represented xinsulin in nM on a log10 scale. The forward simulation of the target system (8) was done by the Euler tau leap method with an integration interval of Δ= 0.001[min]. For the in silico experiment settings, the inter-observation interval was set at Δ=5[min] and T=480 min. That is, we assumed that we could change the concentration of the input insulin every 5 minutes. The observation noise of the in silico experiment was white Gaussian with a variance of 0.0025.

We used xinsulin to indicate the concentration of insulin on a log10 scale and chose Unif [−2, 2] for the uninformative proposal distribution of xinsulin. The lower bound of this distribution (0.01 [nM] = 10−2 [nM]) was based on the basal level of insulin secretion in the work of Kubota et al. [5]. On the other hand, the upper bound (100 [nM]=102 [nM]) was the insulin level at which the response of the target molecules reached saturation (Supplementary Fig. S1). To reflect the constraints, we truncated xinsulin by −2 below and by 2 above. We also assumed that the initial distribution of X(0) was known and set P0 in equation (6) accordingly. We maintained N=1000 particles in each step of the model-based control algorithm. As for the variance parameters, we chose, and . The forward simulation of the discretized ODE mock system was done by the Euler method with an integration interval of Δ= 0.001[min]. We chose t−t−1=5[min]=5000Δ=5000Δ, where Δ and Δ denote the step size of the numerical simulation of the target SDE and our ODE mock system, respectively.

In the control experiment, we repeated the control run of 480 min 50 times, and observed the effects of variation in the initial states across the cell population.

Application to the EGF-dependent Akt pathway in silico

In addition, we applied our model-based control method to a model of the EGF-dependent Akt pathway, developed by Fujita et al. [4] (refer to Supplementary Text S5 for details). In this in silico experiment, we attempted to control the temporal pattern of phosphorylated S6 (pS6) using EGF as an extracellular stimulus (Fig. 2C). That is, stim=EGF and targ= pS6. The forward simulation of the target system was done by the Euler tau leap method with an integration interval of Δ= 0.0001[min]. For the in silico experiment settings, the inter-observation interval was set at Δ=1[min] and T=60 min. That is, we assumed that we could change the concentration of the input EGF every 1 minutes. The remaining settings were the same as in the insulin-dependent AKT pathway.

We used xEGF to indicate the concentration of EGF on a log10 scale and chose Unif [−2, 2] for the uninformative proposal distribution of xEGF. In our in silico experiments, we truncated xEGF by −2 below and by 2 above. As for the variance of the filtering step, we used . The forward simulation of the discretized ODE mock system was done by the Euler method with an integration interval of Δ= 0.0001[min]. We chose t−t−1=1[min]=10000Δ=10000Δ. The remaining settings of the mock system were the same as in the control experiment of the insulin-dependent AKT pathway.

In the control experiment, we repeated the control run of 60 min 50 times, and observed the effects of variation in the initial states across the cell population.

Results

In silico experiments

We applied our model-based control method to the insulin-dependent AKT pathway (Fig. 2A), which was previously modeled by Kubota et al. [5]. Our control objective was to apply an appropriate temporal pattern of the extracellular insulin stimulus to the system in order to make the concentration of a target molecule to trace one of four reference temporal patterns: the pAKT and pGSK3β ramp increase patterns, the G6Pase ramp decrease pattern, and the pS6K triangular pulse pattern. Our choice of the pAKT reference pattern is somewhat special; Kubota et al. could not realize this pattern in their experiments, although they considered it especially important in their simulation study as a pattern that could help investigate the temporal coding in insulin signal transduction.

After some tuning of the future value length τ used in the prediction and control step, our model-based control algorithm was able to make the temporal pattern of intracellular molecules of the in silico system to approximately follow the given reference patterns (Fig. 2B and Supplementary Fig. S2). Similarly, when applied to the EGF-dependent Akt pathway [4] (Fig. 2C), our algorithm was able to make the temporal pattern of the target molecule pS6 closely follow the given ramp increase pattern (Fig. 2D and Supplementary Fig. S3). The results of our in silico experiments suggest that our model-based control method is useful in controlling various intracellular signaling systems, consisting of both multiple parallel pathways (Fig. 2A) and single straight pathway (Fig. 2C).

We examined how the future value length (τ) affects the performance of our model-based control method. We evaluated the control performance in terms of the sum of square errors (SSE) between the reference temporal pattern and the time course of the target intracellular molecule achieved by our model-based control method. A longer future value length improved the control performance for some molecular species, but deteriorated the performance for others (Fig. 3A and 3C). This suggests that the dynamics of the latter group was highly nonlinear and that for this group, long-term predictions were not as reliable as short-term predictions (refer to Discussion for more details).

Figure 3

Control performance. (A) Control performance of the insulin-dependent AKT pathway measured with the sum of square error (SSE) plotted against various choices of the future value length (τ) used in the prediction and control step. Boxplots of SSE between the reference temporal patterns and the trajectories produced in silico by our model-based control method, over 50 different targets (with different initial settings and system noises), are shown. The bottom and top lines of the box represent the first and third quartiles, the black bar the median. The two whiskers flank the 1.5× interquartile range between the upper (third) and lower (first) quartiles. The dots represent outliers. (B) Relationship between the strength of the system noise (σ) and the control performance of the insulin-dependent AKT pathway. (C) Control performance of the EGF-dependent Akt pathway measured with the SSE plotted against various choices of the future value length (τ) used in the prediction and control step. (D) Relationship between the strength of the system noise (σ) and the control performance of the EGF-dependent Akt pathway.

We next examined how the performance of our model-based control method depends on the level of system noise in the target system. We conducted experiments with different values of the parameter σ (0.0, 0.025 (default value), and 0.05) in the SDE of the target system. Because our mock system is noiseless and our method makes no inference for the system noise, the SSE simply increased as the level of the system noise increased (Fig. 3B and 3D). Nonetheless, even with relatively large system noise (σ=0.05), with an appropriate τ setting, the temporal pattern produced by our model-based method did not greatly diverge from the reference pattern (Supplementary Figs. S4 and S5). These results suggest that our model-based control method is robust against the system noise in the target system and hence against the discrepancy between the mock system and the target system.

Comparison with PI controller

Proportional-Integral-Derivative (PID) control is a standard feedback control method. A PID controller monitors the difference between the target time course and the realized time course, and uses the term proportional to the difference (P), the term proportional to the integral of the difference (I), and the term proportional to its derivative (D) in order to determine the control input. Several previous studies applied PI controllers to control problems on biological systems [9,11,14-18,20]. These studies did not use the “D” control, because the level of system and observation noise in their target biological system was large, and “D” controller is known to be sensitive to noise. We therefore compared our model-based control method to a PI controller.

We first determined a set of PI controller parameters that could minimize the SSE between the reference temporal pattern and the temporal pattern produced by the PI controller based on our discretized ODE mock system (Supplementary Text S3 and Supplementary Fig. S6). It should be noted that in our experimental setup the experimenter is only given the mock system and hence, assumed to have only partial knowledge of the target SDE model and the nature of the system noise therein.

When we applied the PI controller to our in silico control problem in the insulin-dependent AKT pathway, in some cases, pAKT and pGSK3β showed oscillatory behaviors (Fig. 4A). On the other hand, pS6K and G6Pase showed stable temporal behaviors in all cases. Stability was also observed for pS6 in the EGF-dependent Akt pathway (Fig. 4C). In terms of SSE, the model-based control method with an appropriate choice of the parameter τ performed significantly better than PI for pAKT (Wilcoxon signed rank test, τ=1, p=1.819*10−5), whereas it performed substantially worse than PI for pS6K (τ=5, p=0.01539) and G6Pase (τ=10, p=2.846*10−8) (Fig. 4B). For pGSK3β (τ=1, p=0.7869) and pS6 in the EGF-dependent Akt pathway (τ=30, p=0.569), there was no statistically significant difference in the performance between our control method and the PI controller (Fig. 4B and 4D). These results suggest that our control method shows superior performance when the prediction of our mock system is accurate, even if its performance is, on average, equal to or not markedly worse than that of the PI controller. Our method also has the advantage of not yielding instable oscillatory behaviors during the control; not a single one of our in silico trials of our algorithm produced such instable behaviors.

Figure 4

Comparison with PI control. (A) Results of the in silico control (σ=0.025) of the insulin-dependent AKT pathway by the PI controller. Red trajectories correspond to the reference temporal patterns, and blue trajectories correspond to the responses of the in silico system. Trajectories of 50 control trials are shown. The variation of the simulated trajectories is due to the system noise and the difference in the initial states. (B) Box-plots of the sum of square error (SSE) between the reference temporal patterns and the controlled trajectories in silico; over 50 different targets are shown. Red, blue, green, and yellow colors correspond to the results of the model-based control method with τ=1, 2, 5 and 10, respectively (‘MB τ=1’, ‘MB τ=2’, ‘MB τ=5’ and ‘MB τ=10’, respectively). The gray color corresponds to the PI control (PI) results. The bottom and top lines of the box represent the first and third quartiles, the black bar the median. The two whiskers flank the 1.5× interquartile range between the upper (third) and lower (first) quartiles. The dots represent outliers. (C) Results of the in silico control (σ=0.025) of the EGF-dependent Akt pathway conducted with the PI controller. (D) Boxplots of the SSE between the reference temporal patterns and the controlled trajectories in silico; over 50 different targets are shown. Red, blue, green and yellow colors represent the results of the model-based control method with τ=1, 10, 20 and 30, respectively (‘MB τ=1’, ‘MB τ=10’, ‘MB τ=20’ and ‘MB τ=30’, respectively). The gray color corresponds to the PI control (PI) results.

Discussion

In this study, we present a nonlinear version of model-predictive control methods with the goal of controlling temporal patterns of intracellular molecules. We verified its performance on the insulin-dependent AKT and EGF-dependent Akt pathways perturbed with system noise. Our model-based control method showed good performance in the in silico experiments. We also confirmed that our algorithm could deal with reference patterns other than the ones presented in Figure 2 (Supplementary Fig. S7). These additional in silico experiments further support the robustness of our methods.

The user should note that the computational cost of our algorithm depends on the dimension and complexity of the mock system that the algorithm uses for prediction and inference, as well as on the future value length τ. For this reason, the effectiveness of our algorithm highly depends on the choice of the mock system. The user should choose a mock system that can mimic the real system as accurately and efficiently as possible.

There are studies that apply the combination between a Monte Carlo filter and model predictive control to non-biological systems, such as an inverted pendulum [24]. In our study, we applied a version of model predictive control to an intracellular signal transduction pathway, a biological system of particular importance in regulating cell behavior. In a similar study, Milias-Argeitis et al. applied a combination between a Monte Carlo filter and model predictive control to a cyanobacteria gene expression system, in which they introduced a mechanism to update online the model parameters based on observations, in order to deal with model uncertainty [19,20]. Here, we conducted model predictive control with a fixed model, but we used the knowledge of the initial distribution of biological molecule concentrations, an aspect of system uncertainty that was not considered in the study by Milias-Argeitis et al. Distributions of intracellular protein concentrations are relatively well studied [25-27].

The dependence of the model-based control method on the choice of the future value length, τ, used in the prediction and control step, is shown in Figure 3. For pAKT and pGSK3β, the control method tended to perform better with a smaller value of τ, whereas for pS6K, G6Pase, and pS6, the control method tended to perform better with a larger value of τ. This suggests that the dynamics of the former group is sensitive to the proximal past and that of the latter group is sensitive to the distant past. This might be related to the speed of signal transduction from the input signal. The insulin or EGF signal is slowly transmitted to S6K, G6Pase, and S6, making their responses sensitive to the distant past. On the other hand, the insulin signal is rapidly transmitted to AKT and GSK3β, making their responses sensitive to the proximal past.

As shown in Figure 4, a portion of cells showed oscillatory behaviors in pAKT and pGSK3β when controlled by the PI controller with conditional integration [31] (see Supplementary Text S3). It might be possible to reduce such oscillatory behaviors by further reducing the inter-observation interval, but this can increase the cost for the prediction and hence is not always a viable option for online control. Our model-based control method achieved stable control with an appropriate choice of the future value length τ. According to a recent study by Fiore et al. [18], the best choice of a control method in terms of performance highly depends on the control task and the complexity of the control target.

In this study, we aimed to control the time course of one target species with the input of one controllable species. In real biological experiments, however, situations may arise in which the researcher needs to control simultaneously the time course of multiple species. This can be extremely difficult, especially when there are nontrivial interactions between the groups of target molecules in the signal transduction pathway. However, when such interactions are weak, and if the series of inputs act on each target molecule “almost” independently, it might be possible to control multiple target species.

The performance of the control method largely depends on the uncertainty of the target biological systems; the system may include system noise, and its initial state may be unknown. Moreover, there might be a notable discrepancy between the target system and the mock system that the algorithm uses for prediction. This is particularly serious for open-loop methods that are not equipped with a mechanism against unexpected behaviors of the system to be controlled. Conversely, the open-loop methods are likely to perform very well in the absence of such uncertainty. For the cellular dynamics studied in [21,22], the open-loop control methods based on ODE were sufficient to achieve satisfactory results. However, some gene expression dynamical systems are known to be highly stochastic. Accordingly, Uhlendorf’s open-loop approach faced difficulties in controlling the average gene expression response of a yeast population to osmotic stress [12]. Milias-Argeitis et al. also reported difficulties in the open-loop control of light-switchable gene expression in yeast [10].

Our in silico experiments assumed that the mock system and the target system coincided on the deterministic/drift component, whereas the nature of the noise was assumed to be unknown. Therefore, the ODE mock system model used the same kinetic parameters as the target model, corresponding to a situation in which the network of the system is reasonably known. However, in real-world biological experiment it might be difficult to obtain the true kinetic parameters of the network.

In order to further confirm the ability of our algorithm to deal with a situation with only partial knowledge of the dynamics, we conducted experiments with a mock system having different kinetic parameters from those of the target system (that is, a case in which our knowledge of both the network and the noise were both incorrect). In particular, we used a mock system with kinetic parameters that were perturbed by the multiplier scale of either 0.9 or 1.1. We showed that, with an appropriate setting of the future value length τ, our algorithm could handle well such a situation (Supplementary Fig. S8).

The variation of parameters and initial states between cells also affects control performance. Carta et al. and Maruthi et al. introduced an algorithm to deal with intrinsic and extrinsic noise of known form [32,33]. They assumed an intrinsic noise that could be expressed in the form of a chemical Langevin approximation and developed a model-based control method of single-cell dynamics of yeast response to osmotic stress. However, it is difficult to identify the form of intrinsic noise in general. In order to derive a chemical Langevin equation that can approximate well the ‘trajectory’ of single-cell dynamics, one needs very specific knowledge of the dynamics at the single-cell level, including the scaling of the reaction rate constants and the relative scale of the number of molecules for each species [34]. Therefore, there is no guarantee that one can improve the control accuracy of single-cell dynamics by using Langevin equations instead of ODEs, unless the detailed form of the system stochasticity is known.

In our approach, we assumed that the experimenter is not given information about the stochasticity and exact initial state of the target biological systems prior to the experiment. In this regard, note that the ODE model we used in our control method agrees only partially with the in silico cellular system that was perturbed by unknown system noise. We would like to emphasize that our in silico model is a nontrivial stochastic system, studied elsewhere [35]. Our in silico experiments were therefore an attempt to control the partially known system, and our success indicates the ability of our algorithm to handle such partially identified systems.

Indeed, unknown overlapping and redundancy of intracellular signal transduction pathways can possibly jeopardize our model-based control method. The extent to which this happens solely depends on the completeness and accuracy of the known mathematical model. The experimental results by Kubota et al. indicate that their mathematical model of the insulin-dependent AKT pathway does not suffer too much from such model indeterminacy [5]. Data driven models [36,37] are also unlikely to suffer much from such indeterminacy. When the destructive effects of model redundancy and overlapping are inevitable, we might be able to improve the control by repeating the online system identification (or, online modification of the model). As discussed above, Milias-Argeitis et al. also adopted a mechanism for changing the parameters of the dynamical system adaptively and updated the model at every observation instance [19,20]. We may also incorporate such a mechanism into our algorithm in order to better handle unexpected behaviors of a system that is only partially known. However, this topic is beyond the scope of the current study.

We might be able to apply our method to in vivo systems. In order to control signal transduction pathways of in vivo systems, we inevitably have to cope with partial knowledge of the target systems. For example, if we are trying to perfectly control in vivo cellular signaling via intravenous stimulus, we first need a perfect model of the dynamics of chemical transport from the vein to individual cells. In general, however, it is impossible to know every detail of the target dynamics, and the dynamics are often perturbed not only by biological noise but also by the noise of our control device. As shown in our in silico experiments, our method is capable of controlling partially known systems. Although the true extent of this capability is yet to be tested, and clearly depends on the target biological system, our work still provides a step toward controlling in vivo systems, for which we are never able to gain complete knowledge. Our work is thus meant to provide an easy-to-implement framework for the control of real biological systems.

Conclusion

In this study, we presented a model-based control method to control the temporal patterns of intracellular signaling molecules using extracellular molecules as stimuli. We validated our method by in silico experiments with unknown system noise. Our method showed good performance in robustly controlling molecule dynamics against unknown system noise of various strengths, even in the absence of complete knowledge of the true model of the target system.

Amongst many dynamical properties of the intracellular signal transduction system, the one that holds particular importance is the relationship between the dynamics of upstream molecules and the dynamics of downstream molecules. A useful tool in understanding this relationship is the methodology to control the dynamics of intracellular molecules with an extracellular stimulus. In this article, we present a model-based control method that dynamically changes the concentration of a controllable upstream molecule to make the target downstream molecule trace a desired time course.