Reduced order modeling and analysis of the human complement system

Abstract

Complement is an important pathway in innate immunity, inflammation, and many disease processes. However, despite its importance, there are few validated mathematical models of complement activation. In this study, we developed an ensemble of experimentally validated reduced order complement models. We combined ordinary differential equations with logical rules to produce a compact yet predictive model of complement activation. The model, which described the lectin and alternative pathways, was an order of magnitude smaller than comparable models in the literature. We estimated an ensemble of model parameters from in vitro dynamic measurements of the C3a and C5a complement proteins. Subsequently, we validated the model on unseen C3a and C5a measurements not used for model training. Despite its small size, the model was surprisingly predictive. Global sensitivity and robustness analysis suggested complement was robust to any single therapeutic intervention. Only the simultaneous knockdown of both C3 and C5 consistently reduced C3a and C5a formation from all pathways. Taken together, we developed a validated mathematical model of complement activation that was computationally inexpensive, and could easily be incorporated into pre-existing or new pharmacokinetic models of immune system function. The model described experimental data, and predicted the need for multiple points of therapeutic intervention to fully disrupt complement activation.

About the Authors:

Adithya Sagar

Roles Software, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, NY, United States of America

Wei Dai

Contributed equally to this work with: Wei Dai, Mason Minot

Roles Investigation, Software, Writing – original draft, Writing – review & editing

Affiliation: Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, NY, United States of America

Mason Minot

Contributed equally to this work with: Wei Dai, Mason Minot

Roles Software, Writing – original draft, Writing – review & editing

Affiliation: Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, NY, United States of America

Rachel LeCover

Roles Software, Visualization

Affiliation: Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, NY, United States of America

Jeffrey D. Varner

Roles Conceptualization, Funding acquisition, Investigation, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, NY, United States of America

ORCID http://orcid.org/0000-0002-2558-7026

Introduction

Complement is an important pathway in innate immunity. It plays a significant role in inflammation, host defense as well as many disease processes. Complement was discovered in the late 1880s where it was found to ‘complement’ the bactericidal activity of natural antibodies [1]. However, research over the past decade has suggested the importance of complement extends beyond innate immunity. For example, complement contributes to tissue homeostasis [2]. It has also has been linked with several diseases including Alzheimers, Parkinson’s, multiple sclerosis, schizophrenia, rheumatoid arthritis and sepsis [3, 4]. Complement also plays positive and negative roles in cancer; attacking tumor cells with altered surface proteins in some cases, while potentially contributing to tumor growth in others [5, 6]. Lastly, several other important biochemical systems are integrated with complement including the coagulation cascade, the autonomous nervous system and inflammation [6]. Thus, complement is important in a variety of beneficial and potentially harmful functions in the body. Despite its importance, there have been few approved complement specific therapeutics, largely because of safety concerns and challenging pharmacokinetic constraints, however, progress is being made [7].

The complement cascade involves many soluble and cell surface proteins, receptors and regulators [8, 9]. The outputs of complement are the Membrane Attack Complex (MAC), and the inflammatory mediator proteins C3a and C5a. The membrane attack complex, generated during the terminal phase of the response, forms transmembrane channels which disrupt the membrane integrity of targeted cells, leading to cell lysis and death. On the other hand, the C3a and C5a proteins act as a bridge between innate and adaptive immunity, and play an important role in regulating inflammation [5]. Complement activation takes places through three pathways: the classical, the lectin and the alternative pathways. The classical pathway is triggered by antibody recognition of foreign antigens or other pathogens. A multimeric protein complex C1 binds antibody-antigen complexes and undergoes a conformational change, leading to an activated form with proteolytic activity. The activated C1-complex cleaves soluble complement proteins C4 and C2 into C4a, C4b, C2a and C2b, respectively. The C4a and C2b fragments bind to form the C4bC2a protease, also known as the classical pathway C3 convertase (C4bC2a). The lectin pathway is initiated through the binding of L-ficolin or mannose binding lectin (MBL) to carbohydrates on the surfaces of bacterial pathogens. These complexes, in combination with mannose-associated serine proteases 1 and 2 (MASP-1/2), also cleave C4 and C2, leading to additional C4bC2a. Thus, the classical and lectin pathways, initiated by different cues on foreign surfaces, converge at C4bC2a. On the other hand, the alternative pathway is activated by a ‘tickover’ mechanism in which complement protein C3 is spontaneously hydrolyzed to form an activated intermediate C3w; C3w recruits factor B and factor D, leading to the formation of C3wBb. C3wBb cleaves C3 into C3a and C3b, where the C3b fragment further recruits additional factor B and factor D to form C3bBb, the alternative C3 convertase (AP C3 convertase) [10]. The role of classical and alternative C3 convertases is varied. First, AP C3 convertases mediate signal amplification. AP C3 convertases cleave C3 into C3a and C3b; the C3b fragment is then free to form additional alternative C3 convertases, thereby forming a positive feedback loop. Next, AP/C4bC2as link complement initiation with the terminal phase of the cascade through the formation of C5 convertases. Both classical and alternative C3 convertases can recruit C3b subunits to form the classical pathway C5 convertase (C4bC2aC3b, CP C5 convertase), and the alternative pathway C5 convertase (C3bBbC3b, AP C5 convertase), respectively. Both C5 convertases cleave C5 into the C5a and C5b fragments. The C5b fragment, along with the complement proteins C6, C7, C8 and multiple C9s, form the membrane attack complex. On the other hand, both C3a and C5a are important inflammatory signals involved in several responses [8, 9]. Thus, the complement cascade attacks invading pathogens, while acting as a beacon for adaptive immunity.

The complement cascade is regulated by plasma and host cell surface proteins which balance host safety with effectiveness. The initiation of the classical pathway via complement protein C1 is controlled by the C1 Inhibitor (C1-Inh); C1-Inh irreversibly binds to and deactivates the active subunits of C1, preventing chronic complement activation [11]. Regulation of upstream processes in the lectin and alternative pathways also occurs through the interaction of the C4 binding protein (C4BP) with C4b, and factor H with C3b [12]. Interestingly, both factor H and C4BP are capable of binding their respective targets while in convertase complexes as well. At the host cell surface, membrane cofactor protein (MCP or CD46) can interact with C4b and C3b, which protects the host cell from complement self-activation [13]. Delay accelerating factor (DAF or CD55) also recognizes and dissociates both C3 and C5 convertases on host cell surfaces [14]. More generally the well known inflammation regulator carboxypeptidase-N has broad activity against the complement proteins C3a, C4a, and C5a, rendering them inactive by cleavage of carboxyl-terminal arginine and lysine residues [15]. Although carboxypeptidase-N does not directly influence complement activation, it silences the important inflammatory signals produced by complement. Lastly, assembly of the MAC complex itself can be inhibited by vitronectin and clusterin in the plasma, and CD59 at the host surface [16, 17]. Thus, there are many points of control which influence complement across the three activation pathways.

Developing quantitative mathematical models of complement could be crucial to fully understanding its role in the body. Traditionally, complement models have been formulated as systems of linear or non-linear ordinary differential equations (ODEs). For example, Hirayama et al., modeled the classical complement pathway as a system of linear ODEs [18], while Korotaevskiy and co-workers modeled the classical, lectin and alternative pathways as a system of non-linear ODEs [19]. More recently, large mechanistic models of sections of complement have also been proposed. For example, Liu et al., analyzed the formation of the classical and lectin C3 convertases, and the regulatory role of C4BP using a system of 45 non-linear ODEs with 85 parameters [20]. Zewde and co-workers constructed a detailed mechanistic model of the alternative pathway which consisted of 107 ODEs and 74 kinetic parameters and delineated between the fluid, host and pathogen surfaces [17]. However, these previous studies involved large models. The central challenge of complement model identification is the estimation of model parameters from potentially limited experimental measurements. Unlike other important cascades, such as coagulation where there are well developed experimental tools and publicly available data sets, the data for complement is relatively sparse. Data sets with missing or incomplete data, and limited dynamic data also make the identification of large mechanistic complement models difficult. Thus, reduced order approaches which describe the biology of complement using a limited number of species and parameters could be important for pharmacokinetic model development, and for our understanding of the varied role of complement in the body.

Materials and methods

Formulation and solution of the complement model equations

We used ordinary differential equations (ODEs) to model the time evolution of complement proteins (xi) in the reduced order model:(1)where denotes the number of reactions and denotes the number of proteins in the model. The quantity τi denotes a time scale parameter for species i which captures unmodeled effects. For the current study, τ scaled with the level of initiator (z) for C5a and C5b; τi = z/z* for i = C5a, C5b where z* was 1 mg/ml, τi = 1 for all other species. The quantity rj(x,ϵ,k) denotes the rate of reaction j. Typically, reaction j is a non-linear function of biochemical and enzyme species abundance, as well as unknown model parameters k (). The quantity σij denotes the stoichiometric coefficient for species i in reaction j. If σij > 0, species i is produced by reaction j. Conversely, if σij < 0, species i is consumed by reaction j, while σij = 0 indicates species i is not connected with reaction j. Species balances were subject to the initial conditions x(to) = xo. Initial conditions were possible were table from Rother [21] and shown in S1 Table.

Rate processes were written as the product of a kinetic term () and a control term (vj) in the complement model. The kinetic term for the formation of C4a, C4b, C2a and C2b, lectin pathway activation, and C3 and C5 convertase activity was given by:(2)where denotes the maximum rate for reaction j, ϵi denotes the abundance of the enzyme catalyzing reaction j, η denotes a cooperativity parameter, and Kjs denotes the saturation constant for species s in reaction j. We used mass action kinetics to model protein-protein binding interactions within the network:(3)where denotes the maximum rate for reaction j, σsj denotes the stoichiometric coefficient for species s in reaction j, and s ∈ mj denotes the set of reactants for reaction j. We assumed all binding interactions were irreversible.

The control terms 0 ≤ vj ≤ 1 depended upon the combination of factors which influenced rate process j. For each rate, we used a rule-based approach to select from competing control factors. If rate j was influenced by 1, …, m factors, we modeled this relationship as where 0 ≤ fij(⋅) ≤ 1 denotes a regulatory transfer function quantifying the influence of factor i on rate j. The function is an integration rule which maps the output of regulatory transfer functions into a control variable. Each regulatory transfer function was modeled using a Hill function. In this study, we used [22]. If a process has no modifying factors, vj = 1. The model equations were implemented in Julia and solved using the CVODE routine of the Sundials package [23, 24].

Model code repository.

The model code, objective functions, experimental data and parameter ensemble used in this study are available under an MIT software license from the Varnerlab GitHub repository [25]. The model equations and kinetic rate expressions are encoded in Balances.jl which is called by the SolveBalances.jl driver function. Note, the user should not directly call SolveBalances.jl. Rather, multiple parameter sets can be simulated by calling the driver function from a script. The kinetic and other model parameters are encoded in DataFile.jl as a dictionary. The parameters stored in this dictionary can be updated in memory to run different simulations. An example script to simulate the model over the parameter ensemble is encoded in sample_ensemble.jl. Plotting routines are encoded in the PlotLib.jl library; these routines plot the experimental data contained in the data subdirectory. The objective functions, which compute the squared residual between the simulations and experimental measurements, are encoded in the complement_lib.jl library. Lastly, we included a set of routines e.g., Make_Fig_2A.jl to solve the model equations, and plot the results versus the experimental data for each simulation figure in this study.

Estimating complement model parameters

A single initial parameter set was estimated using the Dynamic Optimization with Particle Swarms (DOPS) technique [26]. DOPS is a novel hybrid meta-heuristic which combines a multi-swarm particle swarm method with the dynamically dimensioned search approach of Shoemaker and colleagues [27]. DOPS minimized the squared residual between simulated and C3a and C5a measurements with and without zymosan as a single objective. The best fit set estimated by DOPS served as the starting point for multiobjective ensemble generation using Pareto Optimal Ensemble Technique in the Julia programming language (JuPOETs) [28]. JuPOETs is a multiobjective approach which integrates simulated annealing with Pareto optimality to estimate model ensembles on or near the optimal tradeoff surface between competing training objectives. JuPOETs minimized training objectives of the form:(4)subject to the model equations, initial conditions and parameter bounds . The first term in the objective function measured the shape difference between the simulations and measurements. The symbol denotes a scaled experimental observation (from training set j) while the symbol denotes the scaled simulation output (from training set j). The quantity i denotes the sampled time-index and denotes the number of time points for experiment j. The scaled measurement is given by:(5)Under this scaling, the lowest measured concentration become zero while the highest equaled one, where a similar scaling was defined for the simulation output. The second-term in the objective function quantified the absolute error in the estimated concentration scale, where the absolute measured concentration (denoted by ) was compared with the largest simulated value. In this study, we minimized two training objectives, the total C3a and C5a residual w/o zymosan (O1) and the total C3a and C5a residual for 1 mg/ml zymosan (O2). JuPOETs identified an ensemble of N = 2100 parameter sets which were used for model simulations and uncertainty quantification subsequently. JuPOETs is open source, available under an MIT software license. The JuPOETs source code is freely available from the JuPOETs GitHub repository [29].

The simulation and prediction performance of the complement model was measured using the Akaike information criterion (AIC) [30]. In this study, we implemented the AIC as:(6)where Np, Nm denote the number of parameters, and the number of experimental measurements, respectively. The summation term in Eq (6) denotes the residual between the model simulation (x) and experimental measurements (y), where the residual is normalized by the scale of the experimental data (). We compared the AIC for the model parameters estimated in this study, with a random parameter control generated to have a similar order of magnitude. The mean and standard deviation of the AIC was calculated over the parameter ensemble and the random parameter control were reported in this study.

Global sensitivity analysis.

We conducted global sensitivity analysis to estimate which parameters and species controlled the performance of the reduced order model. We computed the total variance-based sensitivity index of each parameter relative to the training residual for the C3a/C5a alternative and C3a/C5a lectin objectives using the Sobol method [31]. The sampling bounds for each parameter were established from the minimum and maximum value for that parameter in the parameter ensemble. We used the sampling method of Saltelli et al. to compute a family of N(2d + 2) parameter sets which obeyed our parameter ranges, where N was the number of trials per parameters, and d was the number of parameters in the model [32]. In our case, N = 400 and d = 28, so the total sensitivity indices were computed using 23,200 model evaluations. The variance-based sensitivity analysis was conducted using the SALib module encoded in the Python programming language [33].

Pairwise sensitivity analysis and clustering.

We perturbed each pair of model parameters by 10% of their nominal value, and then calculated the euclidian distance between the perturbed and nominal system states for physiological conditions. We repeated this calculation for each member of the parameter ensemble, and calculated the mean differences between the perturbed and nominal states. We then clustered the resulting log10 transformed mean distances using the Clustergram routine in MATLAB (The Mathworks, Natick MA). We considered three clusters, high, medium and low displacement.

Robustness analysis.

Robustness coefficients quantify the response of a marker to a structural or operational perturbation to the network architecture. Robustness coefficients were calculated as shown previously [34]. Log-transformed robustness coefficients denoted by were defined as:(7)Here, to and tf denote the initial and final simulation time, while i and j denote the indices for the marker and the perturbation, respectively. A value of , indicates increased marker abundance, while indicates decreased marker abundance following perturbation j. If , perturbation j did not influence the abundance of marker i. In this study, we perturbed the initial condition of C3 or C5 or a combination of C3 and C5 by 50%, 90% and 99% and measured the area under the curve (AUC) of C3a or C5a with and without lectin initiator. We computed the robustness coefficients for a subset of the parameter ensemble (N = 65) and reported the mean robustness value.

Results

In this study, we estimated an ensemble of experimentally validated reduced order complement models using multiobjective optimization. The modeling approach combined ordinary differential equations with logical rules to produce a complement model with a limited number of equations and parameters. The reduced order model, which described the lectin and alternative pathways, consisted of 18 differential equations with 28 parameters. Thus, the model was an order of magnitude smaller than comparable models in the literature. We estimated an ensemble of model parameters from in vitro time series measurements of the C3a and C5a complement proteins. Subsequently, we validated the model on unseen C3a and C5a measurements not used for model training. Despite its size, the model was surprisingly predictive. After validation, we performed global sensitivity and robustness analysis to estimate which parameters and species controlled model performance. Sensitivity analysis suggested CP C3 and C5 convertase parameters were critical, while robustness analyses suggested complement was robust to any single therapeutic intervention; only the knockdown of both C3 and C5 consistently reduced C3a and C5a formation for all cases. Taken together, we developed a reduced order complement model that was computationally inexpensive, and could easily be incorporated into pre-existing or new pharmacokinetic models of immune system function. The model described experimental data, and predicted the need for multiple points of intervention to disrupt complement activation.

Reduced order complement network

The complement model described the initiation of the alternative and lectin pathways, and the downstream integration of the lectin and classical pathways (Fig 1). A trigger event initiated the lectin pathway (encoded as a logical rule), which activated the cleavage of C2 and C4 into C2a, C2b, C4a and C4b, respectively. In this study, we did not explicitly model classical pathway activation mediated by the C1 protein. Instead we described the integration of the lectin and classical pathways at the level of C2 and C4 activation and the formation of the Classical Pathway (CP) C3 convertase (C4aC2b). However, in future studies classical pathway activation could be easily described by simply adding additional terms to the C2/C4 and C2a, C2b, C4a and C4b balances, leaving all other components largely unchanged. Classical pathway (CP) C3 convertase (C4aC2b) then catalyzed the cleavage of C3 into C3a and C3b. The alternative pathway was initiated through the spontaneous hydrolysis of C3 into C3a and C3b. The C3b fragments generated by hydrolysis (or by C4bC2a) could then form the alternative pathway (AP) C3 convertase (C3bBb). However, in this study we considered lumped alternative pathway activation; we did not consider C3w or C3wBb formation. Instead, we formulated a lumped initiation step in which C3a and C3b were formed directly from spontaneous C3 hydrolysis, and factor B was treated as non-limiting. We did this to reduce the complexity both in terms of the number of species as well as the number of parameters that appear in the model. C3b formed from this lumped hydrolysis step then formed the AP C3 convertase. Both the CP and AP C3 convertases catalyzed the cleavage of C3 into C3a and C3b. A second C3b fragment could then bind with either the CP or AP C3 convertase to form the CP or AP C5 convertase (C4bC2aC3b or C3bBbC3b). Both C5 convertases catalyzed the cleavage of C5 into the C5a and C5b fragments. In this study, we simplified the model by assuming both factor B and factor D were in excess. However, we did explicitly account for the action of two other control proteins, factor H and C4BP. Lastly, we did not consider MAC formation, instead we stopped at C5a and C5b. Lectin pathway activation, and C3/C5 convertase activity were modeled using a combination of saturation kinetics and non-linear transfer functions, which resulted in a significant size reduction of the model, while maintaining performance. Binding interactions were modeled using mass-action kinetics, where we assumed all binding was irreversible. Thus, while the reduced order complement model encoded significant biology, it was highly compact consisting of only 18 differential equations and 28 model parameters. Next, we estimated an ensemble of model parameters from time series measurements of the C3a and C5a complement proteins.

[Figure omitted. See PDF.]

Fig 1. Simplified schematic of the human complement system.

The complement cascade is activated through three pathways: the classical, the lectin, and the alternative pathways. Complement initiation results in the formation of classical or alternative C3 convertases, which amplify the initial complement response and signal to the adaptive immune system by cleaving C3 into C3a and C3b. The C3 convertases further react to form C5 convertases which catalyze the cleavage of the C5 complement protein to C5a and C5b. C5b is critical to the formation of the membrane attack complex (MAC), while C5a recruits an adaptive immune response.

https://doi.org/10.1371/journal.pone.0187373.g001

Estimating an ensemble of reduced order complement models

A critical challenge for the development of any dynamic model is the estimation of model parameters. We estimated an ensemble of complement model parameters using in vitro time-series data sets generated with and without zymosan A, a complement pathway activator. The data used for model training was taken from the study of Morad et al. [35] and is given in the Supporting Information (S2 Table). In this study, we assumed zymosan A differentially activated the lectin pathway, consistent with the study of Morad et al [35] and several previous studies [36–38]. However, the role of zymosan may be more complex, as it may activate all three complement pathways under certain conditions. The residual between model simulations and experimental measurements was minimized using the Pareto Optimal Ensemble Technique (JuPOETs), starting from a initial guess generated by the dynamic optimization with particle swarms (DOPS) routine. Unless otherwise specified, all initial conditions were assumed to be at their mean physiological values. While we had significant training data, the parameter estimation problem was underdetermined (we were not able to uniquely determine model parameters). Thus, instead of using the best-fit yet uncertain parameter set, we estimated an ensemble of probable parameter sets to quantify model uncertainty (N = 2100, see Materials and methods). The complement model ensemble captured the behavior of both the alternative and lectin pathways (Fig 2). To estimate alternative pathway model parameters, we used C3a and C5a measurements in the absence of zymosan (Fig 2A and 2B). On the other hand, lectin pathway parameters were estimated from C3a and C5a measurements in the presence of 1 mg/ml zymosan (Fig 2C and 2D). The reduced order model reproduced a panel of alternative and lectin pathway data sets in the neighborhood of physiological factor and inhibitor concentrations. The model fit for parameter sets estimated by JuPOETs, quantified by the Akaike information criterion (AIC), was statistically significantly different than a random parameter control for each case at a 95% confidence level. However, it was unclear whether the reduced order model could predict new data, without updating the model parameters. To address this question, we fixed the model parameters and simulated data sets not used for model training.

[Figure omitted. See PDF.]

Fig 2. Reduced order complement model training.

An ensemble of model parameters was estimated using multiobjective optimization from dynamic C3a and C5a measurements with and without zymosan. The model was trained using C3a and C5a measurements in the absence of zymosan (A–B) or in the presence of 1 mg/ml zymosan (C–D). The solid black lines show the simulated mean value of C3a or C5a for the ensemble, while the dark shaded region denotes the 99% confidence interval of mean. The light shaded region denotes the 99% confidence interval of the simulated C3a and C5a concentration. The experimental training data (points) was taken from Morad et al [35]. All initial conditions not specified by the experimental condition were assumed to be at zero or their physiological serum levels unless otherwise noted (S1 Table).

https://doi.org/10.1371/journal.pone.0187373.g002

We tested the predictive power of the reduced order complement model with data not used during model training (Fig 3). The data used for model validation was taken from the study of Morad et al. [35] and is given in the Supporting Information (S2 Table). Six validation cases were considered, three for C3a and C5a each, respectively. Similar to model training, we compared the AIC for each prediction case to a randomized parameter family. All model parameters and initial conditions were fixed for the validation simulations (with the exception of zymosan, and other experimentally mandated changes). The ensemble of reduced order models predicted the qualitative dynamics of C3a formation (Fig 3, top), and C5a formation (Fig 3, bottom) at three inducer concentrations. For each training case, the AIC was statistically significantly different than the random parameter control for a 95% confidence level. The rate of C3a formation and C3a peak time were directly proportional to initiator dose. Similarly, the C5a plateau and rate of formation were also directly proportional to initiator dose, with the lag time being indirectly proportional to initiator exposure for both C3a and C5a. However, there were shortcomings with model performance. First, while the overall C3a trend was captured (within the 99% confidence interval), the C3a dynamics were too fast with the exception of the low dose case. We believe the C3a time scale was related to our choice of training data, how we modeled the tickover mechanism, and factor B and D limitation. We trained the model using either no or 1 mg/ml zymosan, but predicted cases in a different initiator range; comparing training to prediction, the model performance e.g., the shape of the C3a trajectory was biased towards either high or very low initiator doses. Next, tickover was modeled as a first-order generation processes where C3wBb formation and activity was lumped into the AP C3 convertase. Thus, we skipped an important upstream step which could influence AP C3 convertase formation by attenuating the rate C3 cleavage into C3a and C3b. We also assumed both factor B and factor D were not limiting, thereby artificially accelerating the rate of AP C3 convertase formation. The C5a predictions followed a similar trend as C3a; we captured the long-time C5a behavior but over predicted the time scale of C5 cleavage. However, because the C5a time scale depends strongly upon C3 convertase formation, we can likely correct the C5 issues by fixing the rate of C3 cleavage. Despite these shortcomings, we qualitatively predicted experimental measurements not used for model training typically within the 99% confidence of the ensemble, for three inducer levels. Next, we used global sensitivity and robustness analysis to determine which parameters and species controlled the performance of the complement model.

[Figure omitted. See PDF.]

Fig 3. Reduced order complement model predictions.

Simulations of C3a and C5a generated using 0.1 mg/ml, 0.01 mg/ml, and 0.001 mg/ml zymosan were compared with the corresponding experimental measurements. The solid black lines show the simulated mean value of C3a or C5a for the ensemble, while the dark shaded region denotes the 99% confidence interval of mean. The light shaded region denotes the 99% confidence interval of the simulated C3a and C5a concentration. The experimental validation data (points) was taken from Morad et al [35]. All initial conditions not specified by the experimental condition were assumed to be at zero or their physiological serum levels unless otherwise noted (S1 Table).

https://doi.org/10.1371/journal.pone.0187373.g003

Global analysis of the reduced order complement model

We conducted sensitivity analysis to estimate which parameters controlled the performance of the reduced order complement model. We calculated the total sensitivity of the C3a and C5a residual to changes in model parameters with and without zymosan (Fig 4). In the absence of zymosan (where only the alternative pathway is active), the most sensitive parameter was the rate constant governing the assembly of the AP C3 convertase, as well as the rate constant controlling basal C3b formation via the tickover mechanism. The C5a trajectory was sensitive to the AP C5 convertase kinetic parameters (Fig 4A). Interestingly, neither the rate nor the saturation constant governing AP C3 convertase activity were sensitive in the absence of zymosan. Thus, C3a formation in the alternative pathway was more heavily influenced by the spontaneous hydrolysis of C3, rather than AP C3 convertase activity, in the absence of zymosan. In the presence of zymosan, the C3a residual was controlled by the formation and activity of the C4bC2a, as well as tickover and degradation parameters. On the other hand, the C5a residual was controlled by the formation and activity of CP C5 convertase, and tickover C3b formation in the presence of zymosan (Fig 4B). The lectin initiation parameters were sensitive, but to a lesser extent than CP convertase kinetic parameters and tickover C3b formation. Thus, sensitivity analysis suggested that CP C3/C5 convertase formation and activity dominated in the presence of zymosan, but tickover parameters and AP C5 convertase were more important without initiator. AP C3 convertase assembly was important, but its activity was not. Next, we compared the sensitivity results to current therapeutic approaches; pathways involving sensitive parameters have been targeted for clinical intervention (Fig 4C). In particular, the sensitivity analysis suggested AP/CP C5 convertase inhibitors, or interventions aimed at attenuating C3 or C5 would most strongly influence complement performance. Thus, there was at least a qualitative overlap between sensitivity and the potential of biochemical efficacy. However, total sensitivity coefficients quantify how simultaneous changes in many parameters e.g., rate or saturation constants affect model performance (in this case model fit). To better understand the role of each parameter, and parameter combination, we explored how finite changes in parameter combinations influenced model performance.

[Figure omitted. See PDF.]

Fig 4. Global sensitivity analysis of the reduced order complement model.

Sensitivity analysis was conducted on the two objectives used for model training. A: Sensitivity of the C3a and C5a residual w/o zymosan. B: Sensitivity of the C3a and C5a residual with 1 mg/ml zymosan. The bars denote the mean total sensitivity index for each parameter, while the error bars denote the 95% confidence interval. C: Pathways controlled by the sensitivity parameters. Bold black lines indicate the pathway involves one or more sensitive parameters, while the red lines show current therapeutics targets. Current complement therapeutics were taken from the review of Morgan and Harris [39].

https://doi.org/10.1371/journal.pone.0187373.g004

Pairwise parameter perturbations identified crosstalk within the complement model (Fig 5). We perturbed each pairwise parameter combination by 10%, and calculated the distance between the perturbed and nominal state for each parameter set in the ensemble. We then clustered the mean response of each parameter combination based upon the euclidian distance between the perturbed and nominal states into low (green), medium (red) and high (blue) response clusters. A low response (white) meant the parameter perturbations did not significantly change the system state compared with the nominal case. Four of the 28 parameters (or approximately 14% of the overall model parameters) were in the high response cluster (Fig 5, blue cluster). These parameters included the rate constant controlling the basal formation of C3b (#12), C3a degradation (#26) as well as the catalytic rate constant governing C4bC2a activity (#22). The only C5 related parameter in the high response group was the rate constant controlling the formation of CP C5 convertase (#15). Approximately, 36%, or 10 of the 28 model parameters, were clustered in the medium impact cluster (Fig 5, red cluster). Three parameters (#10, #1, #27) were especially important in this cluster; The reaction order governing C4bC2a activity was important (#10), along with the rate constant controlling C4a and C4b formation from C4 in the lectin initiation pathway (#1), and the constant controlling the inhibitory action of C4BP (#27). Lastly, 50% of the model parameters were clustered in the low response cluster (Fig 5, green cluster). Many of these parameters influenced complement activation; for example, parameter #23 (the C4bC2a saturation constant) was important, just not to the extent of other model parameters. Pairwise synergistic interactions between parameters were also identified. For example, in the high impact cluster, three synergistic relationships were identified, a single positive and two negative cases. Parameters #12 (rate constant governing basal C3b formation) and #15 (formation of CP C5 convertase) acted synergistically to increase the system response. On the other hand, simultaneously changing parameters #12 and #22 or #15 and #26 decreased the system response relative to a single perturbation. However, the most striking examples of synergy occurred in the medium impact cluster; for example, simultaneously increasing parameters #13 (rate constant governing AP C3 convertase formation) and #19 (saturation constant governing AP C5 convertase activity) significantly changed the model state. Changes in parameter #3 (rate constant governing C2a and C2b formation from C2) showed both positive and negative synergistic effects depending upon the other parameter that was perturbed. Taken together, sensitivity coefficients quantified how changes in parameters or parameter combinations affected model performance. However, individual parameters e.g., rate or saturation constants are not easily druggable. To more closely simulate a clinical intervention e.g., administration of anti-complement inhibitors, we performed knock-down analysis on the initial values of C3 and C5 in the absence and presence of flow.

[Figure omitted. See PDF.]

Fig 5. Pairwise sensitivity and clustering of complement model parameters in the presence of 1 mg/ml zymosan.

The response of the complement model was calculated for each parameter combination following a 10% increase in parameter combinations in the presence of 1 mg/ml zymosan. The model parameters were clustered into high (blue), medium (red) and low (green) response clusters based upon the euclidian distance between the perturbed and nominal system state (no perturbation).

https://doi.org/10.1371/journal.pone.0187373.g005

Knock-down analysis in the absence of flow suggested there was no single intervention that inhibited complement activation in the presence of both initiation pathways (Fig 6). Robustness coefficients quantify the response of a protein to a macroscopic structural or operational perturbation to a biochemical network. Here, we computed how the C3a and C5a trajectories responded to a decrease in the initial abundance of C3 and/or C5 with and without lectin initiator. We simulated the addition of different doses of anti-complement inhibitor cocktails by decreasing the initial concentration of C3, C5 or the combination of C3 and C5 by 50%, 90% and 99%. This would be conceptually analogous to the administration of a C3 inhibitor e.g., Compstatin alone or combination with Eculizumab (Fig 4C). The response of the complement model to different knock-down magnitudes was non-linear; a 90% knock-down had an order of magnitude more impact than a 50% knock-down. As expected, a C5 knockdown had no effect on C3a formation for either the alternative (Fig 6A) or lectin pathways (Fig 6B). However, C3a and to a greater extent C5a abundance decreased with decreasing C3 concentration in the alternative pathway (Fig 6A). This agreed with the sensitivity results; changes in AP C3-convertase formation affected the downstream dynamics of C5a formation. Thus, if we only considered the alternative pathway, C3 alone could be a reasonable target, especially given that C5a formation was surprisingly robust to C5 levels in the alternative pathway. Yet, when both pathways were activated, C5a levels were robust to the initial C3 concentration (Fig 6B); even 1% of the nominal C3 was able to generate enough AP/CP C5 convertase to maintain C5a formation. Thus, the only reliable intervention that consistently reduced both C3a and C5a formation for all cases was a knockdown of both C3 and C5. For example, a 90% decrease of both C3 and C5 reduced the formation of C5a by an order of magnitude, while C3a was reduced to a lesser extent (Fig 6B). Taken together, these results suggested that both C3 and C5 inhibitors must be administered in the presence of both initiation pathways.

[Figure omitted. See PDF.]

Fig 6. Robustness analysis of the complement model.

Robustness coefficients were calculated for a 50%, 90% and 99% reduction in C3, C5, or C3 and C5 initial conditions. A: Mean robustness index for C3a and C5a generated in the absence of zymosan. B: Mean robustness index for C3a and C5a generated in the presence of 1 mg/ml zymosan. The color describes the degree of reduction of C3a or C5a following the network perturbation. Robustness coefficients were calculated using all parameter sets with Pareto rank less than five (N = 65). Mean robustness values were reported.

https://doi.org/10.1371/journal.pone.0187373.g006

Discussion

In this study, we estimated an ensemble of experimentally validated reduced order complement models using multiobjective optimization. The modeling approach combined ordinary differential equations with logical rules to produce a complement model with a limited number of equations and parameters. The reduced order model, which described the lectin and alternative pathways, consisted of 18 differential equations with 28 parameters. Thus, the model was an order of magnitude smaller than comparable mathematical models in the literature. We estimated an ensemble of model parameters from in vitro time series measurements of the C3a and C5a complement proteins. Subsequently, we validated the model on unseen C3a and C5a measurements that were not used for model training. Despite its small size, the model was surprisingly predictive. After validation, we performed global sensitivity and robustness analysis to estimate which parameters and species controlled model performance. These analyses suggested complement was robust to any single therapeutic intervention. The only intervention that consistently reduced C3a and C5a formation for all cases was a knockdown of both C3 and C5. Taken together, we developed a reduced order complement model that was computationally inexpensive, and could easily be incorporated into pre-existing or new pharmacokinetic models of immune system function. The model described experimental data, and predicted the need for multiple points of intervention to disrupt complement activation.

There has been a paucity of validated mathematical models of complement pathway activation. To our knowledge, this study is one of the first complement models that combined multiple initiation pathways with experimental validation of important complement products like C5a. However, there have been several theoretical models of components of the cascade in the literature. Liu and co-workers modeled the formation of C3a through the classical pathway using 45 non-linear ODEs [20]. In contrast, in this study we modeled lectin mediated C3a formation using only five ODEs. Though we did not model all the initiation interactions in detail, especially the cross-talk between the lectin and classical pathways, we successfully captured C3a dynamics with respect to different concentrations of lectin initiators. The model also captured the dynamics of C3a and C5a formed from the alternative pathway using only seven ODEs. The reduced order model predictions of C5a were qualitatively similar to the theoretical complement model of Zewde et al., which involved over 100 ODEs [17]. However, we found that the C3a produced in the alternative pathway was nearly three orders of magnitude greater than the C5a generated. While this was in agreement with the experimental data [35], it differed from the theoretical predictions made by Zewde et al., who showed C3a was eight orders of magnitude greater than the C5a concentration [17]. In our model, the time profile of both C3a and C5a generated changed with respect to the quantity of zymosan (the lectin pathway initiator). In particular, the C3a peak time was directly proportional to initiator, while the lag phase for generation was inversely proportional to the initiator concentration. Korotaevskiy et al. showed a similar trend using a theoretical model of complement, albeit for much shorter time scales [19]. Thus, the reduced order complement model performed at least as well as existing larger mechanistic models, despite being significantly smaller.

Global analysis of the complement model suggested potentially important therapeutic targets. Complement malfunctions are implicated in a spectrum of diseases, however the development of complement specific therapeutics has been challenging [3, 39]. Previously, we have shown that mathematical modeling and analysis can be useful tools to estimate therapeutically important mechanisms [40–43]. In this study, we analyzed a validated ensemble of reduced order complement models to better understand the strengths and weaknesses of the cascade. In the presence of an initiator, C3a and C5a formation was sensitive to CP C3/C5 convertase assembly and activity, and to a lesser extent lectin initiation parameters. Formation of the CP convertases can be inhibited by targeting upstream protease complexes like MASP-1,2 from the lectin pathway (or C1r, C1s from classical pathway). For example, Omeros, a protease inhibitor that targets the MASP-2 complex, has been shown to inhibit the formation of downstream convertases [44]. Lampalizumab and Bikaciomab, which target factor B and factor D respectively, or naturally occurring proteins such as cobra venom factor (CVF), an analogue of C3b, could also attenuate AP convertase formation [45–47]. Removing supporting molecules could also destabilize the convertases. Novelmed Therapeutics developed the antibody, NM9401 against propedin, a small protein that stabilizes alternative C3 convertase [48]. Lastly, convertase catalytic activity could be attenuated using small molecule protease inhibitors. All of these approaches are consistent with the results of the sensitivity analysis. On the other hand, robustness analysis suggested C3a and C5a generation could only be significantly attenuated by modulating the free levels of C3 and C5. The most commonly used anti-complement drug Eculizumab, targets the C5 protein [39]. Several other antibodies targeting C5 are also being developed; for example, LFG316 targets C5 in age-related macular degeneration [49], while Mubodina is used to treat atypical hemolytic-uremic syndrome (aHUS) [50]. Other agents such as Coversin [51] or the aptamer Zimura [52] could also be used to knockdown C5. The peptide inhibitor Compstatin and its derivatives are promising approaches for the inhibition of C3 [53]. However, while the knockdown of C3 and C5 affect C3a and C5a levels downstream, the abundance, turnover rate and population variation of these proteins make them difficult targets [54, 55]. For example, the Eculizumab dosage must be significantly adjusted during the course of treatment for aHUS [56]. This suggests that therapeutic targets estimated using whole-body models which incorporate pharmacokinetic factors in combination with biochemistry may give higher fidelity predictions of patient response to therapeutic intervention than static biochemical network models.

Despite its importance, there have been few approved complement specific therapeutics because of safety and pharmacokinetic constraints [7]. A validated complement model, in combination with personalized pharmacokinetic models of immune system function, could be an important development for the field. The integration of effective models of complement and other important cascades in the blood, for example the coagulation cascade, with physiologically based pharmacokinetic models is an exciting opportunity. Physiologically based pharmacokinetic models (PBPK), composed of a collection of organ models (of variable complexity) interconnected by a circulatory system, describe the physical disposition of blood constituents within the body [57]. PBPK models, originally developed to predict the tissue distribution of therapeutic agents or toxins [58], have become an accepted pharmacokinetic tool [59, 60]. PBPK models readily integrate potentially important clinical information such as the demographic and physical characteristics of patients, dissolved oxygen levels, blood flow rates to well- and poorly-perfused organs, pulmonary function and pharmacokinetic factors etc, with blood biochemistry. Thus, PBPK models are ideal candidates to simulate clinically important physical characteristics of patients, while simultaneously simulating biochemical and disease mechanisms on a systems level [61]. We (and many others in the systems pharmacology community) feel this is a critical unmet need, and an opportunity to connect basic science with clinical practice [62].

The performance of the effective complement model was impressive given its limited size. However, there are several questions that should be explored further. A logical progression for this work would be to expand the network to include the classical pathway and the formation of the membrane attack complex (MAC). However, time course measurements of MAC abundance (and MAC formation dynamics) are scarce, making the inclusion of MAC challenging. On the other hand, inclusion of classical pathway activation is straightforward. Liu et al., have shown cross-talk between the activation of the classical and lectin pathways through C reactive proteins (CRP) and L-ficolin (LF) under inflammation conditions [20]. Thus, inclusion of these species, in addition to a lumped activation term for the classical pathway should allow us to capture classical activation. Next, we should address the C3a time scale issue. We believe the C3a time scale was related to our choice of training data, how we modeled the tickover mechanism, and factor B and D limitation. Tickover was modeled as a first-order generation processes where C3wBb formation and activity was lumped into the AP C3 convertase. Thus, we skipped an important step which could strongly influence AP C3 convertase formation by slowing down the rate C3 cleavage into C3a and C3b. The model should be expanded to include the C3wBb intermediate, where C3wBb catalyzes C3 cleavage at a slow rate compared to normal AP or C4bC2as. We also assumed both factor B and factor D were not limiting, thereby artificially accelerating the rate of AP C3 convertase formation. This shortcoming could be addressed by including balances around factor B and D, and including these species in the appropriate kinetic rates. The C5a predictions also had an accelerated time scale. However, because the C5a time scale depended strongly upon C3 convertase formation, we can likely correct the C5 issues by fixing the rate of C3 cleavage. We should also consider including the C2-bypass pathway, which was not included in the model. The C2-bypass mediates lectin pathway activation, without the involvement of MASP-1/2; this pathway could be important for understanding the role of MASP-1/2 inhibitors on complement activation. Lastly, we have assumed that zymosan A differentially actives the lectin pathway, however, this may not be true. We should further explore the potential cross activation of all three branches by zymosan to determine the sensitivity of the model predictions to the assumption of differential activation.

Supporting information

[Figure omitted. See PDF.]

S1 Table. Complement protein initial conditions and associated species indexed used in the model.

All unlisted proteins had an initial condition of zero.

https://doi.org/10.1371/journal.pone.0187373.s001

(PDF)

S2 Table. Experimental measurements used in modeling training and validation from Morad and coworkers [35].

Two sets of C3a and C5a measurements were used for model training while three were used for validation. The information of the experimental data, its usage in this work and its location within the original publication is described below.

https://doi.org/10.1371/journal.pone.0187373.s002

(PDF)

S3 Table. Model species table.

https://doi.org/10.1371/journal.pone.0187373.s003

(PDF)

Acknowledgments

We gratefully acknowledge the suggestions from the anonymous reviewers to improve this manuscript.

Citation: Sagar A, Dai W, Minot M, LeCover R, Varner JD (2017) Reduced order modeling and analysis of the human complement system. PLoS ONE12(11): e0187373. https://doi.org/10.1371/journal.pone.0187373

References

1. Nuttall G. Experimente über die bacterienfeindlichen Einflüsse des thierischen Körpers. Z Hyg Infektionskr. 1888;4:353–394.

2. Ricklin D, Hajishengallis G, Yang K, Lambris JD. Complement: a key system for immune surveillance and homeostasis. Nat Immunol. 2010;11(9):785–797. pmid:20720586

3. Ricklin D, Lambris JD. Complement-targeted therapeutics. Nat Biotechnol. 2007;25(11):1265–75. pmid:17989689

4. Rittirsch D, Flierl MA, Ward PA. Harmful molecular mechanisms in sepsis. Nat Rev Immunol. 2008;8(10):776–787. pmid:18802444

5. Sarma JV, Ward PA. The complement system. Cell Tissue Res. 2011;343(1):227–235. pmid:20838815

6. Ricklin D, Lambris JD. Complement in immune and inflammatory disorders: pathophysiological mechanisms. J Immunol. 2013;190(8):3831–3838. pmid:23564577

7. Ricklin D, Lambris JD. Progress and Trends in Complement Therapeutics. Adv Exp Med Biol. 2013;735:1–22.

8. Walport MJ. Complement. First of two parts. N Engl J Med. 2001;344(14):1058–66. pmid:11287977

9. Walport MJ. Complement. Second of two parts. N Engl J Med. 2001;344(15):1140–4. pmid:11297706

10. Pangburn MK, Müller-Eberhard HJ. The alternative pathway of complement. Springer Semin Immunopathol. 1984;7:163–192. pmid:6238433

11. Walker D, Yasuhara O, Patston P, McGeer E, McGeer P. Complement C1 inhibitor is produced by brain tissue and is cleaved in Alzheimer disease. Brain Res. 1995;675(1):75–82. pmid:7796155

12. Blom AM, Kask L, Dahlbäck B. Structural requirements for the complement regulatory activities of C4BP. J Biol Chem. 2001;276(29):27136–27144. pmid:11369776

13. Riley-Vargas RC, Gill DB, Kemper C, Liszewski MK, Atkinson JP. CD46: expanding beyond complement regulation. Trends Immunol. 2004;25(9):496–503. pmid:15324743

14. Lukacik P, Roversi P, White J, Esser D, Smith G, Billington J, et al. Complement regulation at the molecular level: the structure of decay-accelerating factor. Proc Natl Acad Sci USA. 2004;101(5):1279–1284. pmid:14734808

15. Liszewski MK, Farries TC, Lublin DM, Rooney IA, Atkinson JP. Control of the complement system. Adv Immunol. 1995;61:201–283.

16. Chauhan A, Moore T. Presence of plasma complement regulatory proteins clusterin (Apo J) and vitronectin (S40) on circulating immune complexes (CIC). Clin Exp Immunol. 2006;145(3):398–406. pmid:16907906

17. Zewde N, Gorham RD Jr, Dorado A, Morikis D. Quantitative Modeling of the Alternative Pathway of the Complement System. PloS One. 2016;11(3):e0152337. pmid:27031863

18. Hirayama H, Yoshii K, Ojima H, Kawai N, Gotoh S, Fukuyama Y. Linear systems analysis of activating processes of complement system as a defense mechanism. Biosystems. 1996;39(3):173–185. pmid:8894121

19. Korotaevskiy AA, Hanin LG, Khanin MA. Non-linear dynamics of the complement system activation. Math Biosci. 2009;222(2):127–143. pmid:19854207

20. Liu B, Zhang J, Tan PY, Hsu D, Blom AM, Leong B, et al. A computational and experimental study of the regulatory mechanisms of the complement system. PLoS Comput Biol. 2011;7(1):e1001059. pmid:21283780

21. Rother K, Till G. The complement system. Springer Science & Business Media; 2012.

22. Sagar A, Varner JD. Dynamic Modeling of the Human Coagulation Cascade Using Reduced Order Effective Kinetic Models. Processes. 2015;3(1):178.

23. Bezanson J, Edelman A, Karpinski S, Shah VB. Julia: A Fresh Approach to Numerical Computing. arXiv. 2015;arXiv:1411.1607v4.

24. Hindmarsh AC, Brown PN, Grant KE, Lee SL, Serban R, Shumaker DE, et al. SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers. ACM Trans Math Softw. 2005;31(3):363–396.

25. Varnerlab GitHub Complement Model Repository. https://github.com/varnerlab/Complement_model_repository.git.

26. Sagar A, Shoemaker CA, Varner J. Dynamic Optimization with Particle Swarms (DOPS): A meta- heuristic for parameter estimation in biochemical models. Biotechnol J. 2017; in revision.

27. Tolson BA, Shoemaker CA. Dynamically dimensioned search algorithm for computationally efficient watershed model calibration. Water Res Research. 2007;43(1):W01413.

28. Bassen DM, Vilkhovoy M, Minot M, Butcher JT, Varner JD. JuPOETs: a constrained multiobjective optimization approach to estimate biochemical model ensembles in the Julia programming language. BMC Syst Biol. 2017;11(1):10. pmid:28122561

29. Varnerlab. JuPOETs: A Constrained Multiobjective Optimization Approach to Estimate Biochemical Model Ensembles in the Julia Programming Language;. Available from: https://github.com/varnerlab/POETs.jl [cited 10/06/16].

30. Akaike H. A new look at the statistical model identification. IEEE Trans Auto Cont. 1974;19:716–723.

31. Sobol IM. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comput Simulat. 2001;55:271–280.

32. Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput Phys Commun. 2010;181(2):259–270.

33. Herman, J. SALib: Sensitivity Analysis Library in Python (Numpy). Contains Sobol, Morris, Fractional Factorial and FAST methods. Available online: https://github.com/jdherman/SALib;.

34. Tasseff R, Nayak S, Song SO, Yen A, Varner JD. Modeling and analysis of retinoic acid induced differentiation of uncommitted precursor cells. Integr Biol (Camb). 2011;3(5):578–91.

35. Morad HO, Belete SC, Read T, Shaw AM. Time-course analysis of C3a and C5a quantifies the coupling between the upper and terminal Complement pathways in vitro. J Immunol Methods. 2015;427:13–18. pmid:26391915

36. Rawal N, Rajagopalan R, Salvi VP. Activation of complement component C5: comparison of C5 convertases of the lectin pathway and the classical pathway of complement. J Biol Chem. 2008;283(12):7853–63. pmid:18204047

37. Bergseth G, Ludviksen JK, Kirschfink M, Giclas PC, Nilsson B, Mollnes TE. An international serum standard for application in assays to detect human complement activation products. Mol Immunol. 2013;56(3):232–9. pmid:23787367

38. van Vuuren BJ, Bergseth G, Mollnes TE, Shaw AM. Electroluminescent TCC, C3dg and fB/Bb epitope assays for profiling complement cascade activation in vitro using an activated complement serum calibration standard. J Immunol Methods. 2014;402(1-2):50–6. pmid:24287422

39. Morgan BP, Harris CL. Complement, a target for therapy in inflammatory and degenerative diseases. Nat Rev Drug Discov. 2015;14:857–877. pmid:26493766

40. Luan D, Zai M, Varner JD. Computationally derived points of fragility of a human cascade are consistent with current therapeutic strategies. PLoS Comput Biol. 2007;3(7):e142. pmid:17658944

41. Nayak S, Salim S, Luan D, Zai M, Varner JD. A test of highly optimized tolerance reveals fragile cell-cycle mechanisms are molecular targets in clinical cancer trials. PLoS One. 2008;3(4):e2016. pmid:18431497

42. Tasseff R, Nayak S, Salim S, Kaushik P, Rizvi N, Varner JD. Analysis of the molecular networks in androgen dependent and independent prostate cancer revealed fragile and robust subsystems. PLoS One. 2010;5(1):e8864. pmid:20126616

43. Rice NT, Szlam F, Varner JD, Bernstein PS, Szlam AD, Tanaka KA. Differential Contributions of Intrinsic and Extrinsic Pathways to Thrombin Generation in Adult, Maternal and Cord Plasma Samples. PLoS One. 2016;11(5):e0154127. pmid:27196067

44. Schwaeble HW, Stover CM, Tedford CE, Parent JB, Fujita T. Methods for treating conditions associated with MASP-2 dependent complement activation; 2011. US Patent 7,919,094.

45. Vogel CW, Fritzinger DC, Hew BE, Thorne M, Bammert H. Recombinant cobra venom factor. Molecular immunology. 2004;41:191–199. pmid:15159065

46. Katschke KJ, Wu P, Ganesan R, Kelley RF, Mathieu MA, Hass PE, et al. Inhibiting alternative pathway complement activation by targeting the factor D exosite. J Biol Chem. 2012;287:12886–12892. pmid:22362762

47. Hu X, Holers VM, Thurman JM, Schoeb TR, Ramos TN, Barnum SR. Therapeutic inhibition of the alternative complement pathway attenuates chronic EAE. Mol Immunol. 2013;54:302–308. pmid:23337717

48. Bansal R. Humanized and chimeric anti-properdin antibodies; 2014. US Patent 8,664,362.

49. Roguska M, Splawski I, Diefenbach-Streiber B, Dolan E, Etemad-Gilbertson B, Rondeau JM, et al. Generation and Characterization of LFG316, A Fully-Human Anti-C5 Antibody for the Treatment of Age-Related Macular Degeneration. IOVS. 2014;55:3433–3433.

50. Melis JP, Strumane K, Ruuls SR, Beurskens FJ, Schuurman J, Parren PW. Complement in therapy and disease: Regulating the complement system with antibody-based therapeutics. Mol Immunol. 2015;67:117–130. pmid:25697848

51. Weston-Davies WH, Nunn MA, Pinto FO, Mackie IJ, Richards SJ, Machin SJ, et al. Clinical and immunological characterisation of coversin, a novel small protein inhibitor of complement C5 with potential as a therapeutic agent in PNH and other complement mediated disorders. Blood. 2014;124:4280–4280.

52. Epstein D, Kurz JC. Complement binding aptamers and anti-C5 agents useful in the treatment of ocular disorders; 2007. US Patent App. 12/224,708.

53. Mastellos DC, Yancopoulou D, Kokkinos P, Huber-Lang M, Hajishengallis G, Biglarnia AR, et al. Compstatin: a C3-targeted complement inhibitor reaching its prime for bedside intervention. Eur J Clin Invest. 2015;45(4):423–40. pmid:25678219

54. Sissons J, Liebowitch J, Amos N, Peters D. Metabolism of the fifth component of complement, and its relation to metabolism of the third component, in patients with complement activation. J Clin Invest. 1977;59(4):704. pmid:845257

55. Swaak A, Hannema A, Vogelaar C, Boom F, van Es L, van Aalst R, et al. Determination of the half-life of C3 in patients and its relation to the presence of C3-breakdown products and/or circulating immune complexes. Rheumatol Int. 1982; p. 161–166. pmid:6984977

56. Noris M, Galbusera M, Gastoldi S, Macor P, Banterla F, Bresin E, et al. Dynamics of complement activation in aHUS and how to monitor eculizumab therapy. Blood. 2014; p. 1715–1726. pmid:25037630

57. Kuepfer L, Niederalt C, Wendl T, Schlender JF, Willmann S, Lippert J, et al. Applied Concepts in PBPK Modeling: How to Build a PBPK/PD Model. CPT Pharmacometrics Syst Pharmacol. 2016;5(10):516–531. pmid:27653238

58. Gerlowski LE, Jain RK. Physiologically based pharmacokinetic modeling: principles and applications. J Pharm Sci. 1983;72(10):1103–27. pmid:6358460

59. Sager JE, Yu J, Ragueneau-Majlessi I, Isoherranen N. Physiologically Based Pharmacokinetic (PBPK) Modeling and Simulation Approaches: A Systematic Review of Published Models, Applications, and Model Verification. Drug Metab Dispos. 2015;43(11):1823–37. pmid:26296709

60. Jamei M. Recent Advances in Development and Application of Physiologically-Based Pharmacokinetic (PBPK) Models: a Transition from Academic Curiosity to Regulatory Acceptance. Curr Pharmacol Rep. 2016;2:161–169. pmid:27226953

61. Eissing T, Kuepfer L, Becker C, Block M, Coboeken K, Gaub T, et al. A computational systems biology software platform for multiscale modeling and simulation: integrating whole-body physiology, disease biology, and molecular reaction networks. Front Physiol. 2011;2:4. pmid:21483730

62. Danhof M. Systems pharmacology—Towards the modeling of network interactions. Eur J Pharm Sci. 2016;94:4–14. pmid:27131606