We have collected experimental data from adult C57 mice and built a framework to represent the regulatory relationship among cytokines and macrophages. Based on this framework, we established a set of nonlinear differential equations to characterize the regulatory relationship for macrophage activation in the left ventricle postmyocardial infarction using physical chemistry laws. Our framework and the mathematical model were established based on the following three assumptions.

1)
All monocytes that migrate to the infarct region are differentiated to unactivated macrophages [10].

2)
All activated macrophages are differentiated from unactivated macrophage since previous studies have shown that <5% of macrophages undergo mitotic division [11, 12].

3)
All parameters and coefficients in this model are constant.
Framework of regulatory relationship for macrophage activation
In this framework, myocytes and monocytes were considered as inputs to the system. Cellular densities of M1 and M2 macrophages were considered as the outputs of the system. We chose IL1, IL10, and TNFα as three molecules which regulate macrophage activations in this mathematical model since they were wellrecognized as stimuli for macrophage activation [13–15]. M1 macrophages and myocytes secrete IL1 and TNFα. M2 macrophages secrete IL10 [3]. Further, TNFα and IL1 promote M1 activation and IL_10 promotes M2 activation [16]. IL10 inhibits TNFα, IL1, and itself [17]. The inputoutput and regulating relationship were shown in Figure 1.
Input to the framework
Temporal profiles of monocytes and myocytes densities were used as inputs to our mathematical framework (Figure 1). Myocytes density in healthy adult mice was 6 × 10^{9} cells/ml as an initial value. Myocytes numbers monotonically decreases postMI and is directly associated with LV wall thickness. We have measured the LV wall thickness at days 0, 1, 3, 5, and 7 postMI. The temporal profile of myocytes was determined by combining the initial value and the monotone progression trend (the crosses) as shown in Figure 2(A).
Macrophages density in the left ventricle of healthy adult mice is 2000 cells/ml, which will be used as initial values of unactivated macrophage density in this study [18, 19]. Yang et al has measured temporal profiles of macrophages postMI in mice at days 1, 2, 4, 7, 14, 21, and 28 [20]. The temporal profile of macrophages was obtained by fitting the experimental results as a continuous function as shown in Figure 2(B).
Based on our assumptions 1 and 2, all macrophages were differentiated from monocytes and emigrated from infarct area to the lymph node system. The estimation of unactivated macrophage based on the experimental results [18] is shown as follows,
{M}_{un}\left(t\right)={M}_{un}\left(0\right){e}^{\mu t}+{\int}_{0}^{t}{e}^{\mu \left(ts\right)}Mds,
(1)
where M denotes the differentiation rates of monocytes, M_{
un
} denotes the unactivated macrophage density, and μ denotes the emigration rate of inactivated macrophage. Based on the temporal profile of unactivated macrophage, the monocytes differentiation rate can be estimated as shown in Figure 2(C).
Mathematical model for macrophage activation
The mathematical model of macrophage activations is a set of nonlinear differential equations represented by cellular densities (cell number/ml) of M_{
un
}, M1 and M2, and concentrations (pg/ml) of chemical factors such as IL1, IL10, and TNFα. Cellular densities were determined by the difference between immigration and emigration rates. Concentrations of chemical factors were determined by the balance between their synthesis and degradation rates. The established mathematical model can be written as
\frac{d{M}_{un}}{dt}=\stackrel{differentiation}{\stackrel{\u23de}{M}}\stackrel{activatedtoM1}{\stackrel{\u23de}{{k}_{2}{M}_{un}\frac{I{L}_{1}}{I{L}_{1}+{c}_{IL1}}{k}_{3}{M}_{un}\frac{{T}_{\alpha}}{{T}_{\alpha}+{c}_{{T}_{\alpha}}}}}\stackrel{activatedtoM2}{\stackrel{\u23de}{{k}_{4}{M}_{un}\frac{I{L}_{10}}{I{L}_{10}+{c}_{IL10}}}}\stackrel{emigration}{\stackrel{\u23de}{\mu {M}_{un}}},
(2)
\frac{d{M}_{1}}{dt}=\stackrel{activationeffectsbyI{L}_{1}andTN{F}_{\alpha}}{\stackrel{\u23de}{{k}_{2}{M}_{un}\frac{I{L}_{1}}{I{L}_{1}+{c}_{IL1}}+{k}_{3}{M}_{un}\frac{{T}_{\alpha}}{{T}_{\alpha}+{c}_{{T}_{\alpha}}}}}+\stackrel{transationbetweenM1andM2}{\stackrel{\u23de}{{{k}^{\prime}}_{1}{M}_{2}{k}_{1}{M}_{1}}}\stackrel{emigration}{\stackrel{\u23de}{\mu {M}_{1}}},
(3)
\frac{d{M}_{2}}{dt}=\stackrel{activationeffectsbyI{L}_{10}}{\stackrel{\u23de}{{k}_{4}{M}_{un}\frac{I{L}_{10}}{I{L}_{10}+{c}_{IL10}}}}+\stackrel{transationbetweenM1andM2}{\stackrel{\u23de}{{k}_{1}{M}_{1}{{k}^{\prime}}_{1}{M}_{2}}}\stackrel{emigration}{\stackrel{\u23de}{\mu {M}_{2}}},
(4)
\frac{dI{L}_{10}}{dt}=\stackrel{secretionbyM2}{\stackrel{\u23de}{{k}_{5}{M}_{2}\frac{{c}_{1}}{{c}_{1}+I{L}_{10}}}}\stackrel{degradation}{\stackrel{\u23de}{{d}_{IL10}I\phantom{\rule{0.3em}{0ex}}{L}_{10}}},
(5)
\frac{d{T}_{\alpha}}{dt}=\stackrel{secretionbyM1andmyocytes}{\stackrel{\u23de}{\left({k}_{6}{M}_{1}+\lambda Mc\right)\frac{c}{c+I{L}_{10}}}}\stackrel{degradation}{\stackrel{\u23de}{{d}_{T\alpha}{T}_{\alpha}}},
(6)
\frac{dI{L}_{1}}{dt}=\stackrel{secretionbyM1andmyocytes}{\stackrel{\u23de}{\left({k}_{7}{M}_{1}+\lambda Mc\right)\frac{c}{c+I{L}_{10}}}}\stackrel{degradation}{\stackrel{\u23de}{{d}_{IL1}I{L}_{1}}},
(7)
where M_{
un
}, M_{1}, M_{2} denote the cell densities of unactivated macrophages, M1 macrophages, and M2 macrophages, respectively. Variables IL_{10}, T_{
α
}, and IL_{1} denote the concentrations of IL10, TNFα, and IL1. Variable denotes M the differentiation rate of monocytes and M_{
c
} denotes the myocytes density. The parameters used in these equations with their biological meanings, experimental values, units, and references are listed in Additional file 1. All parameters were determined based on the published data or estimation from other mathematical models [3, 6, 21–28].
Equation 2 determines the density of unactivated macrophages in the infarct area. For the construction part, the unactivated macrophages are differentiated from monocytes as shown in Figure 2(C). For the destruction part, the unactivated macrophages are activated to M_{1} or M_{2}. Additionally, inactivated macrophages do not die locally in the scar tissue but die out in the lymph node system [27].
Equation 3 determines the activation rate of M1 macrophages. For the construction part, IL1 and TNFα promote M1 activation. Parameters k_{2} and k_{3} denote the activation rates of M1 macrophages by IL1 and TNFα. Hill equations are used to represent the promotion effects of IL1 and TNFα and parameters c_{IL 1}and {c}_{{T}_{\alpha}} are the effectiveness of IL1 and TNFα promotion on M1 calculated based on the experimental results [3, 29]. Steinmuller has shown the transition between M1 and M2 phenotypes in vivo [21]. Correspondingly, we use parameter k_{1} to denote the transition from M1 to M2 and parameter {k}_{1}^{\prime} for the transition from M2 to M1 [21]. The destruction part includes emigration of macrophage (μ)and transition from M1 to M2 macrophages (k_{1}) [27].
Equation 4 determines the activation rate of M2 macrophages. The construction part is denoted by activation of M2 macrophages promoted by IL10, and transition from M1 to M2. IL10 promotes M2 activation and this activation rate has been approximated by parameters k_{4} based on the experimental results [3]. The transition rate from M1 to M2 is denoted as k_{1}. The destruction part includes emigration of M2 macrophages (μ) and transition from M2 to M1 macrophages ({k}_{1}^{\prime}), similarly as described in equation 3.
Equation 5 determines the secretion rate of IL10. For construction part, IL10 is secreted by M2 macrophages, and parameter k_{5} denotes the secretion rate of IL10 by M2 macrophages [22, 30]. The destruction part includes the selfinhibition and degradation of IL10. A Hill equation is employed to represent this selfinhibition effect and parameter c_{1} denotes the selfinhibition effect of IL10 postMI [6, 17]. Parameter {d}_{I{L}_{10}} denotes the decay rate of IL10 determined by its halflife time [13].
Equation 6 determines the deposition rate of TNFα, which is secreted by both M1 macrophages and myocytes [6, 23, 28, 31]. We used in vitro results from Meng to determine the secretion rate of TNFα by M1 and secretion rate of TNFα by myocytes (λ) is determined by Horio's experimental results [23, 31]. The inhibition of IL10 is presented by a Hill equation where parameter c represents the effectiveness of IL10 inhibiting TNFα [25]. The destruction part is denoted by the degradation of TNFα. Parameter {d}_{{T}_{\alpha}} is the decay rate of TNFα determined by its halflife time [15, 26].
Equation 7 determines the deposition rate of IL1. IL1 is secreted by both M1 macrophage and myocytes. Parameter k_{7} denotes the secretion rate of IL1 in cultured rat cardiac myocytes [31]. The inhibition of IL10 is presented by a Hill equation similarly as in equation 6 [25]. In the destruction part, parameter d_{IL 1}represents the decay rate of IL1 determined by its halflife time [14].
Stability analysis
If there is no myocardial infarction, monocytes differentiation and myocytes apoptosis should be at a very low level, and the studied macrophage activation pathway should maintain homeostasis. We have calculated the equilibrium point of the system without any input and performed Lyapunov stability analysis. Our analysis showed that without any monocytes differentiation and myocytes secretion, the system would stay at the origin when t → ∞.
In the case of myocardial infarction, myocytes apoptosis and necrosis triggered inflammatory responses and significant monocytes differentiation, which will drive the system to a new equilibrium point. Correspondingly, the cell densities of M1 and M2 increase postMI. We have obtained a steady state as E = [20, 1200, 3500, 0.73, 1.1, 5.9] from our computational simulations. The steady states match with the experimental measurements collected from healthy adults without myocardial infarction [32]. In addition, the stable region of the established mathematical model depends on the strength of the input.
Computational results
Computational simulations of macrophage activation were carried out by solving the nonlinear differential equations with MATLAB. The initial conditions of unactivated, M1 and M2 macrophage densities were chosen as M_{
un
} = 2000 cells/ml and M_{1}(0) = M_{2}(0) = 0 cells/ml. The concentrations of IL1, TNFα, and IL10 were set as 0.1 pg/ml. The inputs of this system were shown in Figure 2. Outputs of the system, M1 and M2 densities, were shown in Figure 3. The concentrations of IL1, IL10 and TNFα were shown in Figure 4. Our computational results were shown as solid lines while the experimental results were shown as discrete crosses in these figures [17, 33, 34].
Our computational results demonstrated that from days 0 to day 3 postMI, cellular densities of the M1 phenotype increased at a faster rate than the M2 phenotype. At day 10, the M2 phenotype dominates over the M1 phenotype. This prediction agrees with the results reported by Troidl [4]. Additionally, temporal profiles of IL1 and TNFα significantly increased from days 0 to 1 postMI in our computational simulations, which match the experimental results reported by Sumitra [33]. Comparison between the computational and experimental results demonstrates a similar trend for the temporal profiles, suggesting the effectiveness of our model.