SENSITIVITY OF SIGNALING PATHWAY DYNAMICS TO PLASMID TRANSFECTION AND ITS CONSEQUENCES

. This paper deals with development of signaling pathways models and using plasmid-based experiments to support parameter estimation. We show that if cells transfected with plasmids are used in experiments, the models should include additional components that describe explicitly eﬀects induced by plasmids. Otherwise, when the model is used to analyze responses of wild type, i.e. non-transfected cells, it may not capture their dynamics properly or even lead to false conclusions. In order to illustrate this, an original mathematical model of miRNA-mediated control of gene expression in the NF κ B pathway is presented. The paper shows what artifacts might appear due to experimental procedures and how to develop the models in order to avoid pursuing these artifacts instead of real kinetics.


1.
Introduction. Responses to any external stimuli, behavior and fate of any living cell is determined by interplay of many complex biochemical processes that take place inside it. One of the most important of them is gene transcription, in which information stored in DNA is transcribed into newly produced messenger RNA (mRNA) molecules, that subsequently are used as templates to produce new proteins. These proteins, in turn, are the active players in all intracellular processes.
The process of gene transcription is regulated by one or more mechanisms, involving various molecules. Among them, so called micro RNA (miRNA) has recently gained a lot of interest. Those short non-coding RNA molecules may bind to specific sequences in mRNA, which results in either mRNA degradation or translation blocking. There are many species of miRNA, each being able to target various mRNA species.
Knowledge of intracellular regulatory mechanisms helps to advance our understanding of the nature of many illnesses [2] as well as provides hints for finding better therapies or new therapeutic targets. However, revealing only the structure of these mechanisms is not sufficient, as the dynamics of the processes involved plays an important role in cellular responses to changes in environment, stress factors or DNA damage. Therefore, mathematical modeling of regulatory pathways that control intracellular biological and chemical processes is a rapidly expanding field in biomedical research [19,11,4], supported by developments in new experimental techniques (e.g., [3]).
2.1. miRNA-based control mechanisms. There are many mechanisms involved in regulation of gene expression. Some of them utilize small, noncoding RNA molecules called micro-RNAs (miRNAs). These molecules are implicated in many biological processes [20,7] and are the subject of ongoing research in cancer and other diseases [10]. The human genome encodes over a thousand miRNA species [5,13], targeting most protein coding genes [8,12].
Processes leading to activation of miRNAs and detailed mechanisms of their actions are very complex and describing them is beyond the scope of this paper. Here, we consider only the final result of their activity, which in most cases is the reduction of the number of proteins coded by the mRNA that contains a specific Binding Site (BS) for a given miRNA. This is achieved either by targeting the mRNA for degradation or by prohibiting ribosomes to attach to mRNA and start the translation process. Gaining knowledge of the sequences of the BS facilitated development of experimental procedures, in which artificial plasmids containing these specific BSs are introduced into cells, making it possible to observe the results of miRNA control actions.

2.2.
Experimental procedures involving transfection. In many cases experimental research on miRNA is based on so called plasmid transfection [9,17]. In such procedures, a reporter gene containing a sequence of interest and coding a bioluminescent reporter (specific for the miRNA that is a subject of the investigation) is introduced into an expression vector, which is subsequently transferred into cells. Following transfer, the cells are assayed for the presence of the reporter by directly measuring the reporter protein itself or the enzymatic activity of the reporter protein [1]. The idea is to measure the activity without disturbing other processes that take place inside the cell.
Plasmids used in the experiments under consideration contain specific BSs, which are the same as in the target gene of interest, controlled by a specific miRNA. As a result, the miRNAs of interest will bind to mRNA coded by the plasmid, as well as to the mRNA that they are supposed to control. Observation of luminescence of plasmid-coded proteins allows to measure miRNA activity, as binding of miRNA to plasmid-coded mRNA leads to reduction of the number of bioluminescent reporter proteins.
3. A mathematical model of a miRNA-controlled regulatory modulethe case study.

Consequences of transfection experiments.
Let us first consider a generic system, in which a specific miRNA regulates expression of mRNA of interest. Binding of miRNA to its BS in the mRNA results in creation of a complex, which after some time dissociates releasing miRNA and driving mRNA to degradation ( Fig. 1 with solid line elements only). In a wild type cell, in the simplest case (no feedback between mRNA and miRNA in their production), these processes may be described by the following set of three ODEs: where m 1 , r 1 and c 1 denote concentrations of mRNA, free (unbounded) regulatory miRNA and mRNA-miRNA complex, respectively. Parameters p m , p r , d m , d r denote the rates of mRNA production, mRNA degradation, miRNA production and miRNA degradation, respectively, while k 1 and k −1 are complex association and dissociation rates, respectively.
If the experiments used to measure results of miRNA activity involve plasmid transfection, the system under consideration contains also other components, as depicted in Fig. 1. Then, the equation (2) should include additional terms, describing processes showed in Fig 1 in dashed lines, thus changing into where m 2 and c 2 are concentrations of plasmid-coded mRNA containing BSs for miRNA and its complex with regulatory miRNA, respectively. Their dynamics is governed by the following equations: If the association and dissociation rate constants k 1 and k −1 are estimated from plasmid based experimental data, the estimation procedure must be based on the full model that includes (4), (5) and (6). Moreover, transfection is not a fully controlled process in the sense that the number of plasmid copies entering each cell is a random variable. As a result, different cells obtain a different number of plasmid copies. There are two implications of this fact. First, the production rate p m2 in (5) may be different in each cell. Second, it is the average luminescence that is measured in the experiments, so to compare simulation and experimental results one should first run a large number of simulations with p m2 being the random parameter, and then average these results before comparing them to experimental data.
3.2. The model of miRNA-based regulation of NFκB pathway. So far, most of the models describing control mechanisms of intracellular processes have been focused on activation or inhibition of gene transcription. Relatively few works deal with modeling of miRNA-mediated regulation. Among them, interesting works can be found, devoted to analysis of either small regulatory modules or larger interaction networks (e.g. [21,23,24,25]).
However, none of them takes into account the nature of experiments whose results the model output should reflect. The simple generic example presented in the previous section indicates that this may lead to at least poor parameter estimation. In order to check if that is the case in modeling a specific signaling pathway, we have built an original model that includes miRNA-mediated regulation in the NFκB pathway. The model is based on interaction between two subsystems that can be found in literature. The first of them is one of the simplest regulatory modules, in which levels of intracellular species, including NFκB are controlled by specific miRNA molecules [23], as more complex models (e.g. [24]) require introduction of molecular species that cannot be measured experimentally.
Its graphical representation is shown in Fig. 2, in the right dashed line-bounded box. It consists of two interacting factors, NFκB and IL-6 that induce the production of each other, as well as two miRNAs. Each of these factors also induces production of a specific miRNA, that inhibits the expression of the other factor.
It should be noted that the original equations from [23] had to be modified to provide nonnegativness of solutions for any set of parameters. In [23] the effect exerted by miRNA on NFκB was given by a linear term. While for the parameters chosen in that paper it appeared that such model provides an excellent fit to experimental data, the solutions are highly dependent on the absolute level of miRNA which is significantly altered in transfection experiments. For sufficiently large miRNA values, such simplification lead to negative concentrations of NFκB.
The second subsystem, depicted in Fig. 2 in the left dashed line-bounded area is a well established model of NFκB pathway [14]. It has been chosen as the most representative one. Though there exist other models of this pathway (e.g. [6,18,22]), they are all based on the same feedback loops. Since this paper aims at showing the approach to model development and not the analysis of a particular pathway, Figure 2. Block diagram of the miRNA-based regulation of the NFκB pathway. Two subsystems, described in the literature (adapted from [14,23]) are distinguished by the dashed lines. The molecule p21 represents plasmid coding mRNA that contains miRNA-specific binding sites, introduced in the transfection experiment. An arrow between two elements indicates that increase in one molecular species leads to increase of the other. A line section at the end of the graph edge represents an inhibitory action.
an arbitrary choice of the pathway model does not constrain its applicability in any sense.
The complete model is described by a set of ordinary differential equations, that describes the changes in molar concentrations of molecular players involved in the pathway, with the variables named after the molecules (Table 1). All modifications to the original equations from [23] and [14] have been introduced in the standard way, following the approach based on the law of mass action.
The parameters appearing in these equations are explained in Table 2. The parameter original values have been modified, as adding new interactions is required to obtain the same system responses when no plasmids are added. It should be noted that the model represents only one mode of miRNA actions, consisting in transcript degradation. If it was aimed at capturing regulatory mechanism of translation inhibition, the equation should include another term, representing dissociation of miRNA-mRNA complexes. Nevertheless, even when it is the degradation action that is modeled, these complexes should be also included explicitly as variables, which will be proved in the Results section.
Unfortunately, there exist no publicly available data showing changes in concentrations of main molecular players both in transfected and untransfected cells.
free cytoplasmic IκBα protein (IκBα) n free nuclear IκBα protein plazmid-coded mRNA with miR21 binding site (P 21 ) protein coded by the plasmid mRNA Table 1. Model variables Therefore, to illustrate the differences between models built on distinct assumptions, only simulation results are shown.

4.
Simulation results and discussion.
4.1. Simulation procedure. Transfection efficiency vary from one cell to another in the sense that each cell contains a different number of plasmid copies. This can be incorporated in the model by randomizing parameter n in (23) and repeating the simulation for a predetermined number of times. Each simulation represents behavior of one cell, whose intracellular processes are governed by deterministic equations with only one initial condition being a random variable. Sample simulation results are shown in Fig. 3, in which four distinct cases are presented. The first plot illustrates the behavior of a wild type cell, which has not been transfected with plasmids. Basically, this is the system that is of an interest and whose dynamics is to be reflected by the model. Three other cases represent the behavior of cells in a population that is transfected with the plasmids: with high efficiency (n has been sampled from the uniform distribution on the [20,30] interval), low efficiency (n has been sampled from the uniform distribution on the [5,10] interval) or variable efficiency (in which some cells do not contain the plasmid at all, other might contain many copies -n has been sampled from the uniform distribution on the [0, 30] interval). The uniform distribution has been chosen arbitrarily as the most relevant in this case. Since n should be a nonnegative integer, and in a realistic situation about 30 copies is the maximum number that can be accepted by a cell. Thus, using other distributions, more commonly used in systems biology, seems not justified. IκB-NFκB association rate k 7 inducible IκBα mRNA production rate c 5 miR21-NFκB association rate k 8 miR21 inducible production rate c 10 miR146-IL association rate k 9 inducible miR146 production rate d 1 IKKn and IKKi degradation rates k 10 inducible IL production rate d 2 IκBα degradation rate k 11 inducible plasmid mRNA production rate d 3 NFκB degradation rate p 1 neutral IKK production rate d 4 A20 mRNA degradation rate p 2 constitutive IKK production rate d 5 A20 protein degradation rate p 3 constitutive NFκB production rate d 6 IκBα protein degradation rate p 4 constitutive A20 mRNA production rate d 7 IκBα mRNA degradation rate p 7 constitutive IκBα mRNA production rate d 8 miR21 degradation rate p 8 constitutive miR21 production rate d 9 miR146 degradation rate p 9 constitutive miR146 production rate d 10 IL degradation rate p 10 constitutive IL mRNA production rate d 11 plasmid mRNA degradation rate p 11 constitutive plasmid mRNA production rate d 12 plasmid degradation rate t 5 A20 translation rate e 6 IκBα nuclear export rate t 6 IκBα translation rate e 7 IκBα-NFκB nuclear export rate t 12 plasmid translation rate k v cytoplasmic to nuclear volume ratio Table 2. Parameters used in simulation.
To compare simulation results with experimental data, one must consider the experimental setup again. In addition to luminescence measurements of the plasmidcoded protein, Reverse Transcription Polymerase Chain Reaction (RT-PCR) technique and Western Blots are used to measure transcript and protein content of a   given molecular species, respectively. RT-PCR provides information about the fold change in mRNA copies while Western Blotting produces images, in which light intensity of stripes in specific locations is proportional to the protein amount. These techniques are averaging over cell populations. Moreover, in most cases it is the average luminescence that is used as an indicator of how much of the reporter protein is present in cells. This of course leads to artificial attenuation of the oscillation, which is not present at the single cell levels but is nevertheless a widely used approach in many experiments. Therefore, in his paper we assumed such procedures are used in experimental data acquisition. The obvious result of the introduction of many plasmid copies into a cell is an increased competition for miRNA between its original targets and plasmid-coded mRNA. Due to this competition, a smaller number of original targets is bound to miRNA and more of them escape miRNA-based regulation. Consequently, the more plasmids are in a cell, the less regulation is observed at the level of miRNA targets (see NFκB levels in Fig. 4). However, specific normalization procedures (like normalization by the area under the curve -AUC) are required when comparing simulation and experimental data [15]. Therefore, such normalization has been included in the calculations whose results are shown in this paper (Fig. 4b) and, consequently, it is the qualitative behavior that is the subject of the discussion. 4.2. Change of system dynamics by plasmid transfection. Introducing plasmids to the cell in transfection experiments changes the dynamics of system responses to external stimuli, not only in quantitative but also in a qualitative way. There is shift in oscillation phase visible in the NFκB time course (Fig. 4(b)). While the initial system response looks the same regardless of transfection, the differences become visible in later stages. If the measurement procedure is adjusted so that even smaller oscillations are captured, and both model and measurements are normalized to facilitate their comparison, the influence of transfection efficiency is apparent. The oscillation phase is shifted, with two consequences: (i) without the plasmid variable, the model would be fitted to specific experiment conditions which do not appear in wild type cells; (ii) if the time series consists of only several time points (e.g. 0.5, 1, 1.5, 2 hrs), the dynamics of mRNA regulated by miRNA, indicated by them may be completely different.
What is even more important is that the phase of oscillations of plasmid-coded mRNA may be completely different depending on transfection efficiency, as illustrated in (Fig. 3). Time course of some protein level may indicate increasing trends, which is due only to the effects introduced by the experimental procedure. It does not make sense to try to fit a model of wild type cell behavior to such trends. Instead, it can be retraced after the model with plasmids are fitted to data from experiments with transfection and, subsequently, the plasmid levels has been set to zero (Fig. 5), i.e. initial condition for (p21 t ) variable and the parameters n and c 11 .  It is clear from Figs. 4-5 that the qualitative dynamics is different in transfected cells and normal cells. Both the oscillation frequencies and their damping ratios differ from one cell type to another. If a model was built without taking into account the experimental procedure, its behavior would not reflect the dynamics of cells that normally do not contain artificial plasmids.
Another interesting observation concerns the oscillation amplitude in system responses. Numerical experiments suggest that they depend on the transfection efficiency (Fig. 6). Proteins that are coded by miRNA-regulated transcripts may exhibit much larger oscillations in transfected cells. Time course of miRNA concentration in wild type cells, shown in black line in Fig. 6 can be retraced after the model with plasmids are fitted to data from experiments with transfection and, subsequently, the plasmid levels has been set to zero. Otherwise, the model for wild type cells would exhibit oscillations on a scale not existing in these cells. 4.3. miRNA-mRNA complexes as necessary variables. Before fitting numerical results to experimental data, careful review of experimental procedure is required to recognize the actual signal that is measured. Often, the RT-PCR measurements are reported as a direct source of valid biological conclusions. However, in RT-PCR procedures it is the total amount of mRNA that is measured, including both molecules accessible for ribosome complexes and those bound to miRNA and hence inactivated. Only the former is the template for protein production. So, the effects of miRNA regulatory actions may not be visible in RT-PCR experiments. Then, experimental results may prove difficult to interpret with almost constant mRNA and changing protein levels. However, knowing the intricacies of experimental procedures, one can modify the model to take them into account.
In the case the miR146-mediated regulation of IL-6 was under consideration, an additional variable should be introduced, representing concentration of miR146 complexes with IL-6, whose dynamics is governed by the following equation: Then, the system output that should be compared to the results of RT-PCR experiments, is given by the total mRNA present in the system, i.e. the sum of free mRNA and miR146-IL complexes. It is the mathematical model then, that provides an added value to analysis of experimental results, making it possible to discern free mRNA serving as a template for protein production. The conclusions should therefore be drawn after retrieving its fraction from the total amount (Fig. 7). Fig. 7 shows that when the regulatory mechanisms lead to oscillations in mRNA levels, they may not be visible in RT-PCR experiments. Moreover, if the measurement procedure is adjusted so that even smaller oscillations are captured, and both model and measurements are normalized to facilitate their comparison, the influence of transfection efficiency on the dynamics of the observed variable becomes apparent. In the first one, creation of mRNA-miRNA complexes is followed by their degradation. In the second case, the complexes are simply used to inhibit translation, and eventually the mRNA may return to the original, unbound state. In most cases, it is not known which mode of action is utilized in specific miRNA-mRNA systems. Appropriate changes in the mathematical model may help in answering that question. Te model presented in the preceding sections was based on the assumption that miRNA-mRNA complexes are degraded following their creation. If a simple inhibition was the true mode of action, equations (21) and (22) should be modified, to include terms describing dissociation of miR146 complexes with IL-6: Equation (25) describing changes in concentration of the complex should be left unchanged, but the interpretation of the d 13 parameter would be different. It would be a dissociation rate constant, instead of the degradation rate constant.
Following these changes, behavior of two alternative models can be compared. While the levels of miRNA and plasmid protein might be exactly the same in both models (if the value of dissociation and degradation rate constants, represented by d 13 in both models was the same), the total plasmid mRNA allows to distinguish these cases. While this conclusion might have been reached without mathematical modeling, simulation allows to explain why plasmid mRNA and protein oscillations are observed to be in phase in some experiments and completely off-phase in others. Once again, it is the varying level of plasmid transfection that is responsible for differences between experiments (Fig. 8) (d) Figure 8. Influence of the transfection efficiency on mRNA measurements: (a) total plasmid mRNA when miRNA action leads to degradation; (b) total plasmid mRNA when miRNA action leads to inhibition without degradation; (c) plasmid protein level; (d) free miRNA level. Red, green and blue lines represent low, high and variable transfection efficiency. Black line in (d) represents miRNA level in wild type cells 4.5. Implications for experiment planning. If measurements taken at discrete time points are used to fit the model to data, the model kinetics would be quite different from the actual one. In transfected cells the NFκB level increases between 2 and 4 hours, whereas if the plasmid is turned off in the model, representing wildtype cell behavior, it decreases in that period of time (Fig. 9).  Figure 9. Dynamics exhibited by discrete measurements of miRNA-regulated mRNA. Black, red, green and blue points represent system response in wild-type cells, transfected with low, high and variable efficiency, respectively.
The differences in the oscillation frequencies seem to be the most important. One of the practical applications of mathematical modeling in analysis of intracellular processes is in experiment design. Due to high costs of experiments, amount of molecular species is usually measured only at several time points (unless live cell imaging is the experimental technique). These are often chosen arbitrarily, most often as the same time points that can be found in literature, to facilitate better comparison with previous experiments. However, if the nature of experiments is different (i.e., experiments with and without transfection) the time points should be chosen differently. In Fig 9 one can see that the maximum and minimum values of the concentrations in transfected and non-transfected cells are off-phase. Of course, it is possible to fit parameters of a model that does not include transfection elements so that the dynamics it exhibits reflects the observations. However, if one attempts to subsequently work with such model to analyze the behavior of the wild-type cells, which is the main goal of research in this area, it might prove to be useless. For example, if the aim was to control the cell fate by targeting molecules at the peak of their activity, the resulting control action would be exerted at precisely the moment of their lowest activity.

5.
Conclusions. In this paper we have shown that a mathematical model describing dynamics of intracellular processes should be developed with particular attention to experimental procedures used to obtain the data to which the model is fitted. Numerical analysis presented in the paper indicates that if experimental procedures involve cells transfection with plasmids, these plasmids should be explicitly included in the model. Moreover, such models should include the concentration of miRNA-mRNA complex as a separate variable whose dynamics is traced. Experiments cannot focus on mRNA measurement only, since conclusions drawn from them may be misleading. This is particularly important in the short time-horizon experiments, as the protein levels may exhibit oscillations that are off-phase comparing to mRNA oscillations.
Another important observation stemming from this paper is that proteins that are coded by miRNA-regulated transcripts may exhibit much larger oscillations in transfected cells. Therefore, when conclusions drawn from transfection experiments are transposed on behavior of wild type cells, careful analysis of simulation results is needed.
An additional contribution of this paper is a novel model of miRNA-based regulation of the NFκB pathway. However, it should not be regarded as the most important one, as it requires further experimental validation.