Algorithmic multiscale analysis for the FcRn mediated regulation of antibody PK in human.

A demonstration is provided on how algorithmic asymptotic analysis of multi-scale pharmacokinetics (PK) systems can provide (1) system level understanding and (2) predictions on the response of the model when parameters vary. Being algorithmic, this type of analysis is not hindered by the size or complexity of the model and requires no input from the investigator. The algorithm identifies the constraints that are generated by the fast part of the model and the components of the slow part of the model that drive the system within these constraints. The demonstration is based on a typical monoclonal antibody PK model. It is shown that the findings produced by the traditional methodologies, which require significant input by the investigator, can be produced algorithmically and more accurately. Moreover, additional insights are provided by the algorithm, which cannot be obtained by the traditional methodologies; notably, the dual influence of certain reactions depending on whether their fast or slow component dominates. The analysis reveals that the importance of physiological processes in determining the systemic exposure of monoclonal antibodies (mAb) varies with time. The analysis also confirms that the rate of mAb uptake by the cells, the binding affinity of mAb to neonatal Fc receptor (FcRn), and the intracellular degradation rate of mAb are the most sensitive parameters in determining systemic exposure of mAbs. The algorithmic framework for analysis introduced and the resulting novel insights can be used to engineer antibodies with desired PK properties.

Introduction

Pharmacokinetics (PK) modeling is an indispensable tool in the drug development process, where mathematical models are implemented to characterize and predict the drug distribution and concentration over time profile[1]. Owning to the advance in computational and analytical methods, as well as the increasing need for preclinical-to-clinical translation and dosing optimization, the PK models are transforming from simple one and two compartmental models to more mechanistically based models. These models usually consist of a large number of complex equations and a much larger number of parameters, making it difficult to assess the impact of specific parameters, reactions and variables on a certain output[2-4]. The multiscale character of these models introduces another difficulty in assessing the influence of the various parameters and reactions, since the available in vivo data are usually collected in clinical trials with a frequency that adheres to the slow time scales[5]. As a result, processes that are characterized by fast time scales (such as the antibody-receptor binding) are not adequately resolved, so that identifiability problems arise[6]. Model reduction has been recognized as a viable remedy to the obstacles created by the complexity of the models of interest and their multiscale character. Several methodologies are available for model reduction, such as sensitivity analysis, singular perturbation (i.e., time-scale analysis), lumping, optimization etc; see Refs.[2,7-11] for details.

The model reduction methods employed most frequently are the quasi-steady-state approximation (QSSA) and partial-equilibrium approximation (PEA); e.g.,[6,12-18]. Both methodologies are based on the multiscale character of the full model and are limiting cases of singular perturbation analysis[19-22]. The major drawback in their use is that their validity requires the identification of the number of fast time scales, the fast variables and fast reactions. Such identifications depend strongly on the value of the parameters employed and on the period of interest along a specific trajectory. These drawbacks are all manifested in the well known Michaelis-Menten (MM) model, for which the sQSSA, rQSSA and tQSSA reduced models were developed, each valid in distinct domains of the parameter space[23-25]. Depending on the parameters, the validity of these reduced models might be restricted in only a portion of a specific trajectory, so that more than one models must be employed if it is desired to approximate the full trajectory[26]. The drawbacks of QSSA and PEA and the fact that their use is not free of ambiguities are all well known. As a result, their implementation must always be subjected to verification[27-29].

The traditional paper-and-pencil singular perturbation analysis (SPA) faces similar problems QSSA and PEA do; i.e., the number of fast variables and the proper non-dimensional form of the full model must be identified, both of which are strongly dependent on the value of the parameters employed or on the period of interest along a specific trajectory. However, a number of algorithms were developed that can produce most or all SPA does. All these algorithms are based on the geometric singular perturbation analysis (GSPA), which was developed in an effort to systematize the various SPA methodologies[30-32]. According to GSPA, the fast dynamics of a multiscale model are responsible for the development of low-dimensional surfaces in phase space, on which the process is constrained to evolve according to the slow dynamics[33]. The algorithm that provides all SPA does is the Computational Singular Perturbation (CSP), which identifies (1) the number of fast time scales and the reactions that are responsible for their development, (2) the variables that can be considered fast, (3) the constraints that develop due to fast dynamics and the reactions that participate in the related equilibria and (4) the reactions that drive the process within the constraints[34-36]. CSP is used extensively in combustion and chemical kinetic studies and recently in problems in biology and pharmacokinetics[26,37-43]. Regarding the subject of the present investigation, CSP has been already employed for the analysis of two target-mediated drug disposition (TMDD) models; a one-compartment model and a two-compartment with first order absorption one[40,41,43]. Among others, it was shown that algorithmic CSP tools provide (1) system-level understanding and (2) all necessary information to predict the response of the system when a change in the parameters is introduced.

Here, a simple antibody model, considering Neonatal Fc Receptor (FcRn) recycling of Immunoglobulin G (IgG), will be analyzed with CSP in order to demonstrate how system-level understanding can be obtained algorithmically. IgG typically shows an extended half-life compared to other protein molecules of similar size (i.e., 150 kDa) due to neonatal Fc receptor (FcRn) mediated recycling[44,45]. The binding between FcRn and IgG is negligible at physiological pH but strong at the acidic pH within the endosomes. This pH-dependent binding protects IgG from degradation in the lysosomal space, while allowing to release it back into the system circulation upon return to the cell surface[46]. Modification of the affinity between FcRn and IgG can alter antibody serum half-life as desired[47-49]. Based on the importance of FcRn binding on IgG half-life, mechanistic-based PK models of IgG usually incorporates FcRn binding component. The IgG enters the endosomal space via pinocytosis, and the binding of FcRn in endosomal space can salvage IgG from degradation. While this process has been described in the past, a systemic analysis of each parameters involved in this process has not been performed before. Consequently, here we have employed the CSP algorithm to gain more insight into the FcRn mediated recycling process of antibodies. The present investigation demonstrates the ability of CSP to provide insight that cannot be provided by the conventional methods.

The structure of the manuscript is as follows. First, the antibody model will be introduced and two QSSA-based reduced models will be generated, which are the appropriated ones for the parameter set employed here and are frequently encountered in the related literature[47,48,50,51]. A methodology for improving the accuracy of these two QSSA models will be presented and the physical understanding provided by the QSSA and the more accurate models will be compared. The CSP formalism will be briefly stated and the CSP analysis will follow. It will be demonstrated how CSP diagnostics can (1) provide system-level understanding and (2) suggest ways to modify the parameters of the model according to a desired response of the system.

Model

A schematic representation of the model employed here is depicted in Fig. 1. This model consists of two compartments: a plasma compartment and endosomal compartment. In the model, IgG can move from plasma compartment to endosomal compartment via pinocytosis (CL). IgG can then interact with FcRn via association () and dissociation () rate constants. Bound IgG will recycle back to plasma space (CL) and unbound IgG in the endosomal space will be eliminated with a elimination rate constant (). The model was able to capture the wildtype and Fc-mutated IgG PK profile in human after IV administration; e.g.,[46]. The process is governed by the system:whereand and . and denote IgG concentration in the plasma and the endosom, respectively, denotes the FcRn concentration in the endosom and denotes the IgG.FcRn concentration in the endosom. From the last two relations of Eq. (1), it follows that the conservation law holds. The system in Eq. (1) can be cast in the compact form:where and denote the stoichiometric vector and rate of the i-th reaction, respectively.

Figure 1

Schematic representation of the model formulating the interaction of IgG with the FcRn receptor, , , , ; .

Conventional reduced models

The conventional methodology that is frequently employed for the construction of reduced models will be employed next. In particular, the QSSA for and will be considered, which have previously employed in similar investigations; e.g.,[47,48,50,51]. These assumptions are the ones appropriate for the parameter set that was employed for the numerical demonstrations that follow. Of course, different assumptions might be appropriate if different parameter sets are considered.

Reduced models based on the QSSA of and

Suppose that reaction 3f: is very fast and that its reactant IgG is the fast variable, so that the following QSSA for holds:Substituting in Eq. (1), by eliminating the rate of the fast reaction 3f, yields:A reduced model can be constructed, using the algebraic relation in Eq. (5) and the differential Eqs. (6, 8, 9). Note that Eqs. (8, 9) keep satisfying the conservation law , which can replace either Eq. (8) or (9).

Suppose now that the 4th reaction is very fast and that its reactant IgG.FcRn is the fast variable, so that the following QSSA for holds:Substituting in Eqs. (6–9), by eliminating the rate of the fast reaction 4, yields:A new reduced model can be constructed, using the algebraic relations in Eqs. (5) and (10), the differential Eq. (11) and the conservation law .

Higher order correction to the reduced models based on the QSSA of and

A higher order correction for the QSSA of is obtained by differentiating Eq. (5) with time[52,53], substituting from Eq. (1) and solving for the rate of the fast reaction :orwhere and . Substituting in Eq. (1) yields:A reduced model can be constructed by using Eqs. (17–20) or by using the algebraic relation Eq. (15) (which must be solved for the fast variable ) and the differential Eqs. (17, 19, 20). Note that this reduced model simplifies to the one in Eqs. (6–9) in the limits and .

Similarly, a higher order correction for the QSSA of is obtained by differentiating Eq. (10) with time[52,53], substituting from Eqs. (17–20) and solving for the rate of the fast reaction :orwhere , and then substituting in Eqs. (17–20):A reduced model can now be constructed by using Eqs. (23–26) or by using the algebraic relations Eqs. (15) and (21) (which must be solved for the fast variables and ) and the differential Eqs. (23, 25). Note that this reduced model simplifies to the one in Eqs. (11–14) in the limits , and .

The conservation law is still valid and can replace either of Eqs. (19,20) and either of Eqs. (25,26).

Discussion

At this point it must be emphasized that for the construction of the reduced models in Eqs. (5–9), (10–14), (15–20) and (21–26) the investigator must identify (1) the proper number of QSSAs and (2) the variables and reactions that are considered fast. These are identifications, which are usually based on the experience of the investigator and must be subject to an a posteriori verification.

A comparison of (1) the reduced model in Eqs. (5–9), generated with one QSSA, with the equivalent one in Eqs. (15–20) and of (2) the reduced model in Eqs. (10–14), generated with two QSSAs, with the equivalent one in Eqs. (21–26), reveals that the models in Eqs. (15–20) and (21–26) provide more informations than those in Eqs. (5–9) and (10–14), respectively. For example, the differential equation for in Eq. (7) does not indicate the way the various reactions influence the evolution of this concentration. In contrast, the equivalent Eq. (18), indicates that the evolution of is determined by the 1st, 2nd and 4th reactions. This influence of the 4th reaction is absent in the related to equation of the full model, Eq. (1). However, it is included in Eq. (18) because its influence is exercised via the equilubrium in Eq. (15). Similar comments apply for the equations governing the evolution of and .

It will be shown next that the reduced models in Eqs. (15–20) and (21–26) are more accurate than those in Eqs. (5–9) and (10–14). Given a complex model, the analytical derivation of such more accurate reduced models requires significant input from the investigator and considerable mathematic manipulations. However, as it will be shown next, this type of reduced models can be constructed algorithmically with the CSP algorithm. In addition, it will be demonstrated that the CSP-generated models provide a full system-level understanding, in contrast to the models in Eqs. (15–20) and (21–26). Such an understanding, is most useful when considering changes in the model or when exploring different operating conditions for the function of the model.

CSP formalism and diagnostics

In CSP context, the 4-dim. model in Eq. (4) is cast in the form:where the state vector is defined as , the 4-dim. column vectors are the CSP basis vectors and are the related amplitudes:where the 4-dim. row vectors are the CSP dual basis vectors (), and and are the stoichiometric vector and rate of the i-th reaction introduced in Eq. (4). The term in Eq. (27) represents a CSP mode and is characterized by a single time scale, say . The CSP modes are ordered in Eq. (27) according to the value of ; first the fastest mode (smallest ), etc. Due to the conservation law , and .

CSP generated reduced models

When the fastest mode becomes exhausted, Eq. (27) reduces to:Similarly, when the two fastest modes become exhausted, Eq. (27) reduces to:In general, considering an N-dim. system, when M fast modes become exhausted, the reduced model is:The M relations define the constraints (i.e., the low dimensional surface in phase space) on which the system evolves according to the reduced model that is expressed by the differential equation in Eq. (31).

CSP diagnostic tools

CSP provides a number of algorithmic tools that allow system level understanding. A detailed description of the CSP tools is provided in the Supplement. In summary, the following tools are available: Given a solution from a model in the form of Eq. (1), the CSP diagnostics (API, II, TPI and Po) can be computed algorithmically on the basis of the expressions provided in the Supplement, with relatively simple codes in Fortran, Python, Mathematica, MatLab, etc. The Fortran code employed here is publicly available in GitHub; see Additional Information. The required input are (1) the number of variables N and of unidirectional reactions K, (2) the stoichiometric vectors , the reaction rates and their gradient .

Amplitude Participation Index (API) Identifies the reactions that contribute significantly to the formation of each of the M constraints, expressed by the relations () in Eq. (31). The related index is , which measures the contribution of the k-th reaction to the cancellations occurring in the expression . is scaled so that , where K is the number of reactions in the model[35,37,54].

Importance Index (II) Identifies the reactions contributing the most in driving the evolution of each variable according to the system of differential equations in Eq. (31). The related index is , which measures the influence of the k-th reaction in driving the n-th variable, within the established constraints. is scaled so that [35,37]. Positive (negative) indicates an influence of the k-th reaction in increasing (decreasing) the n-th variable.

Time scale Participation Index (TPI) Identifies the reactions that contribute significantly to the development of the time scale that characterises the n-th mode (). The related index is measures the contribution of the k-th reaction to the n-th time scale and is scales so that [37,55,56]. A negative (positive) implies that the k-th reaction contributes to a dissipative (explosive) character of the n-th timescale . By definition, dissipative (explosive) timescales relate to the components of the system that tend to drive it towards (away from) equilibrium[34,35].

Pointer (Po) Identifies the variables related the most to the generation of the n-th mode. The related index is and is scaled so that [35,54,57]. A relatively large value of indicates that the i-th species is strongly associated to n-th CSP mode and the n-th timescale. A value of close to unity suggests the validity of the Quasi Steady-State Approximation (QSSA) for the i-th variable, while only two large values of indicate the validity of the Partial Equilibrium Approximation (PEA)[21]. The variables that are pointed by the Pointer of the M exhausted modes are those that are “slaved” to the remaining variables via the equilibrium relations () in Eq. (31) and are those that adjust the most when the related equilibrium relation is perturbed[54].

CSP results

The model parameters are presented in Table 1. Physiological values of human plasma and endosomal volume are used for the two compartments separately. The association () and dissociation () rate constant between IgG and FcRn and the IgG endosomal degradation rate () were obtained from the literature[49,58]. The pinocytosis and exocytosis rate () was calculated based on the reported per-unit endosomal uptake rate multiplying the endosomal volume.

Table 1

The parameter set employed for modelling the system in Eq. (1).

V_p=3.10L	k_on=0.559(nM·h)^-1	k_off=23.9h^-1	k_deg=25.0h^-1
V_e=0.34L	CL_up=0.167L·h^-1	FcRn_0=4.98×10^4nM	Dose=24.0×10^4nmol

Considering the parameters in Table 1, the solution of the model in Eq. (1) was analyzed with CSP and the results are displayed in Figs. 2, 3, 4 and 5. Figures 2 and 3 show that the profiles of (1) the concentrations, (2) the reaction rates and (3) the concentrations rate of change, exhibit periods of rapid changes and periods of very small variations. It will be shown that the latter periods are indicative of established equilibria.

Figure 2

The evolution of the variables , , and (left) and of the rates (k = 1, 2, 3f, 3b, 4) (right). The periods indicated and refer to those in which one and two modes, respectively, are exhausted; i.e., in the period (Eq. (29)) and in the period and (Eq. (30)).

Figure 3

The evolution of the rate of change of , , and .

Figure 4

The evolution of the time scales (left) and of the amplitudes (right); .

Figure 5

Evolution of CSP diagnostics for the three modes, ; TOP: Pointers ; , , , . MIDDLE: Time scale Participation Indices ; . BOTTOM: Amplitude Participation Indices ; .

Figure 4 displays the evolution of the timescales and of the amplitudes of the first three modes in Eq. (27) (); the forth mode is neglected, because it represents the conservation law , for which and . It is shown that significant gaps exits between and and between and , so this feature established the multiscale character of the model. At each point in time the characteristic time scale is the fastest of those having a non-negligible amplitude; i.e., in the period the characteristic time scale is (M = 0 in Eq. (31)) and no reduced model is valid, in the period the characteristic time scale is (M=1) and the reduced model in Eq. (29) is valid and in the period the characteristic time scale is (M=2) and the reduced model in Eq. (30) is valid.

In Figs. 2, 3, 4 and 5 the periods in which a reduced model is valid is denoted by an arrow; indicates the period in which one mode is exhausted () and a model of the form of Eq. (29) is valid, while indicates the period in which two modes are exhausted ( and ) and a model of the form of Eq. (30) is valid.

The M = 1 period

Consider the period, in which the first equilibrium holds. According to Fig. 4, the fastest time scale that characterizes the development of this equilibrium there is much faster than ; about a four orders of magnitude difference. It is also shown that, right before the start of this period, the amplitude decreases considerably, as the first mode becomes exhausted, while the amplitude of the second mode remains at very large values. Figure 5 shows that the variable related to the first mode in the period is () and that the reaction responsible for the emergence of is the IgG()-consuming reaction 3f in the endosom (). Finally, it is shown that the equilibrium expressed by the relation involves the equilibration of , first with (; and ), then with () and finally with (; and ). During the period, the dominant mode is the second one, the time scale of which is mainly generated by reaction 4 and by reaction 2 at the last part of the period, as shown in Fig. 5. Indicative values for all CSP diagnostics for this period () are provided in Table 2, as they were computed at .

Table 2

Solution data and CSP diagnostics at , where .

y_i (nM)	CpIgG = 77377.76	CeIgG = 1.70	CeFcRn = 49422.84	CeIgG.FcRn = 377.16
RHS (nMh^-1)	gpIgG = − 4148.09	geIgG = 33.84	geFcRn = − 37744.56	geIgG.FcRn = 37744.56
τ_i (h^-1)	τ_1=0.36 10^-4	τ_2=1.77	τ_3=458.09
f^i	f^1 = 1.86 10^-3	f^2 = 3.78 10^4	f^1 = 1.53 10^2
n	IgG,p	IgG,e	FcRn,e	IgG.FcRn,e
Dn1	0.000	0.999	0.000	0.001
Dn2	0.091	0.000	0.001	0.906
Dn3	0.908	0.000	0.036	0.055
k	1	2	3f	3b	4
R^k (nMh^-1)	38006.135	42.479	46943.920	9014.109	185.251
Jk1	0.000	- 0.001	- 0.998	- 0.001	0.000
Jk2	- 0.086	- 0.034	0.032	- 0.032	- 0.817
Jk3	- 0.199	- 0.211	0.203	- 0.199	0.187
ck1	0.99910488	- 0.99910293	- 1.00000080	1.00000080	0.00089592
ck2	1.00209268	- 0.90652448	0.00080178	- 0.00080178	- 1.00289446
ck3	0.00403610	0.09556626	- 0.00008647	0.00008647	- 0.00394963
Pk1	0.404	- 0.000	- 0.500	0.096	0.000
Pk2	0.992	- 0.001	0.001	- 0.000	- 0.005
Pk3	0.941	0.025	- 0.025	0.005	- 0.004
	r_1 = 1.94 10^-6	r_2 = 8.98 10^-4	a = 0.1096

All these informations lead to the conclusion that in the period the QSSA of the variable is valid; i.e., . In order to assess the accuracy of such an approximation, the following relative errors were computed:where is the relative error in approximating generated by the QSSA expression in Eq. (5), is the relative error generated by the higher order correction in Eq. (15) and is the relative error generated by the CSP derived expression:where, . The results displayed in Table 3 indicate that the expression in Eq. (15) provided higher order accuracy to the standard QSSA in Eq. (5) () and the CSP-derived Eq. (35) provided higher to Eq. (15) accuracy (). Note that using the values displayed in Table 2, the expression for the the first equilibrium, as derived by CSP in Eq. (35) at , yields:which is similar to the QSSA relation in Eq. (5) and its extension in Eq. (15). The values of coefficients in this expression stay practically constant throughout the period. In addition, according to the values displayed in Fig. 5 , the major participants in the equilibrium across this period are , and ; i.e., and provide negligible contributions and can be neglected. Therefore, Eq. (35) can be simplified as:providing good accuracy across the period.

Table 3

Relative error in approximating at indicative points along the period: via the QSSA expression in Eq. (5) (), the more accurate expression in Eq. (15) () and the CSP expression in Eq. (35) ().

	t (h)	10^-2	10^-1	1	2	3	4
R^3f-R^3b	er11	0.89×10^-3	0.98×10^-3	0.29×10^-2	0.11×10^-1	0.19×10^-2	0.86×10^-3
	er21	0.39×10^-5	0.42×10^-5	0.98×10^-5	0.33×10^-4	0.58×10^-4	0.62×10^-4
	er31	0.50×10^-7	0.58×10^-7	0.47×10^-6	0.56×10^-5	0.11×10^-6	0.21×10^-7

The M = 2 period

Consider now the period, in which the two equilibria and hold. Regarding the first equilibrium, the CSP diagnostics are similar to those in the period; (1) the pointed variable is reactant of the fast reaction 3f, which is the main responsible for the emergence of and (2) the fact that and that both indicate that the QSSA for is still valid, as in the period.

Regarding the second equilibrium Fig. 4 shows that the time scale that characterizes its development is much faster than ; about two orders of magnitude. It is also shown that the amplitude decreases considerably right before the start of the period, as the second mode becomes exhausted, while the amplitude of the third mode remains at very large values. Figure 5 shows that at the start of this period the variable related the most to this mode is , followed by , while for the remaining part of the period is only the latter variable (). Figure 5 also shows that the reaction responsible for the emergence of is initially the IgG()-consuming reaction 2 and then the IgG-FcRn()-consuming reaction 4. Finally, Fig. 5 shows that the equilibrium expressed by the relation involves mainly the equilibration of and with (), in which contributes only at the start of the period. Indicative values of all CSP diagnostics for this period are provided in Table 4, as they were computed at .

Table 4

Solution data and CSP diagnostics at , where .

y_i (nM)	CpIgG = 31346.03	CeIgG = 58.84	CeFcRn = 21202.58	CeIgG.FcRn = 28597.42
RHS (nMh^-1)	gpIgG = − 148.07	geIgG =− 0.58	geFcRn = 120.40	geIgG.FcRn = − 120.40
τ_i (h^-1)	τ_1=0.84 10^-4	τ_2 = 1.53	τ_3 = 100.82
f^i	f^1 = 0.55 10^-6	f^2 = 0.37	f^1 = 148.03
n	IgG,p	IgG,e	FcRn,e	IgG.FcRn,e
Dn1	0.000	0.995	0.003	0.002
Dn2	0.068	0.004	0.096	0.831
Dn3	0.932	0.000	0.476	− 0.407
k	1	2	3f	3b	4
R^k (nMh^-1)	15369.432	1471.039	697404.378	683478.404	14046.381
Jk1	0.000	−  0.002	−  0.996	0.000	−  0.002
Jk2	−  0.060	−  0.149	0.060	−  0.061	−  0.668
Jk3	−  0.290	−  0.258	0.109	−  0.107	0.234
ck1	0.99525512	−  0.99525062	−  1.00000977	1.00000977	0.00475466
ck2	1.00855828	−  0.92519772	0.00190065	−  0.00190065	−  1.01045892
ck3k	0.01881054	0.08335607	−  0.00017575	0.00017575	−  0.01863478
Pk1	0.011	−  0.001	−  0.499	0.489	0.000
Pk2	0.461	−  0.040	0.039	−  0.038	−  0.421
Pk3	0.316	0.134	−  0.134	0.131	−  0.286
	r_1 = 4.53 10^-6	r_2 = 4.78 10^-3	a = 0.1096

All these informations lead to the conclusion that in the period the QSSA of the and variables are valid. In order to assess the accuracy of these approximations, the relative errors , and in approximating () in Eqs. (32–34) were computed along with the relative errors:where is the relative error in approximating generated by the QSSA expression in Eq. (10), is the relative error generated by the higher order correction in Eq. (21) and is the relative error generated by the CSP derived expression:where, . The results displayed in Table 5 indicate that the expressions in Eqs. (15) and (21) provided higher order accuracy to the standard QSSA in Eq. (5) and (10) ( and ) and the CSP-derived Eqs. (35) and (41) provided similar or better to Eqs. (15) and (21) accuracy ( and ).

Table 5

Relative error in approximating (top) and (bottom) at indicative points along the period: TOP: via the QSSA expression in Eq. (5) (), the more accurate expression in Eq. (15) () and the CSP expression in Eq. (35) (), BOTTOM: via the QSSA expression in Eq. (10) (), the more accurate expression in Eq. (21) () and the CSP expression in Eq. (41) ().

	t (h)	5	10	20	50	100	1000
R^3f-R^3b	er11	0.92×10^-3	0.70×10^-3	0.42×10^-3	0.12×10^-3	0.41×10^-4	0.44×10^-5
	er21	0.59×10^-4	0.47×10^-4	0.32×10^-4	0.16×10^-4	0.99×10^-5	0.42×10^-5
	er31	0.24×10^-7	0.15×10^-7	0.55×10^-8	0.43×10^-9	0.39×10^-10	0.24×10^-13
R^4	er12	0.88×10^-2	0.97×10^-2	0.10×10^-1	0.10×10^-1	0.86×10^-2	0.44×10^-2
	er22	0.87×10^-4	0.48×10^-4	0.13×10^-4	0.12×10^-3	0.10×10^-3	0.18×10^-4
	er32	0.22×10^-4	0.23×10^-4	0.28×10^-4	0.12×10^-3	0.26×10^-4	0.24×10^-7

As Table 5 indicates, the accuracy in approximating with the CSP-derived expression is higher than that in approximating . This feature relates to the fact that the accuracy of the constraint , which is generated after the M-th time scale becomes exhausted (so that the characteristic time scale changes from to ), is proportional to the size of the ratio ; i.e., the smaller this ratio the higher the accuracy[59,60]. Accordingly, the accuracy of the CSP expression for approximating in Eq. (35), which becomes valid after the fastest time scale becomes exhausted, is proportional to , while the accuracy of the CSP expression for approximating in Eq. (41), which becomes valid after becomes exhausted, is proportional to . Figure 4 shows that . Therefore, the accuracy of the former approximation will be higher than the latter, in accordance to the results in Table 5. Note that using the values displayed in Table 4, in the expressions for the first and second equilibria, as derived by CSP in Eqs. (35) and (41) at , yield:As in the period discussed previously, Eq. (42) is still similar to the QSSA relation in Eq. (5) and its extension in Eq. (15). Equation (43) is similar to the QSSA relation in Eq. (10) and resembles to its extension in Eq. (21), in which is however not present. The values of coefficients in these two expressions vary negligibly in the M = 1 period. In addition, according to the values displayed in Fig. 5 (1) the major participants in the first equilibrium () across the period are and , while , and provide negligible contributions and (2) the major participants in the second equilibrium () are and , while , and provide small but no negligible contributions only at the early stage of this period. As a result, the previous expressions can be simplified throughout the period as:where the terms in the brackets provide smaller contributions.

The reactions driving the process

When the M equilibria have been established, as it is expressed by the algebraic relations (), the system is driven by the dominant reactions in the simplified system of differential equations in Eq. (31). These reactions are identified by the Importance Index, introduced in “CSP diagnostic tools” Section and displayed in Fig. 6; () implies that the variable is produced (consumed) by the reaction.

Figure 6

The participation index () for all four variables (,,,) in the periods where .

Figure 6 shows that initially (while ), is consumed only by reaction 1 (), is produced by reaction 1 and then gradually consumed by reaction 3f () and and are gradually consumed and produced, respectively, by reaction 3f. These findings can also be predicted by inspection of the governing equations in Eq. (1).

After the first equilibrium is established (: ) and the 1st mode becomes exhausted, the process is driven by the remaining 2nd and 3rd modes; see Eq. (29), which is valid in the period. Figure 6 shows that the action of the two active modes in this period refers mainly to reactions 1 () and 4 () and partly to reaction 2 (), while the contributions of reactions 3f and 3b manifest towards the end of this period cancelling each other. In particular, Fig. 6 shows that is consumed by reaction 1 () and is produced by reaction 4 (). and are produced by reaction 1 (, ) and during the last stage of this period are consumed mainly by reaction 4 (, ) and partly by reaction 2 (, ). These reactions have exactly the opposite influence to the evolution of (mainly , and partly ). At the end of this period, reactions 3f and 3b exhibit a minor influence to the evolution of the three variables in the endosom, but their contributions cancel each other.

When the second equilibrium is established, both the 1st and 2nd modes are exhausted (: ) and it is only the 3rd mode that drives the system; see Eq. (30), which is valid in the period. It is shown in Fig. 6 that in this period reactions 1 and 2 tend to decrease , and (, , , , , ) and to increase (, ), reaction 4 has the opposite tendencies (, , , ), while the contributions of reactions 3f and 3b cancel each other.

Some of the features displayed in Fig. 6 can be detected from the original model in Eq. (1) in the period where no mode is exhausted (), the reduced models in Eqs. (6–9) or Eqs. (17–20) that are valid in the period and the reduced models in Eqs. (11–14) or (23–26) that are valid in the period; e.g., the influence of reactions 1, 3f and 4 during the initial and the periods or the transition from reaction 3f to reaction 1 as the major driving reaction for the evolution of and at the start of the period.

However, there are features that cannot be foreseen by inspection of the full or the reduced models; even in the case where these models are simple, as those considered here. For example, Fig. 6 shows that during the transition from the to the period, the influence of reactions 1 and 4 reverses. That is, in the period reaction 1 tends to increase and and to decrease (, , ), while in the period its influence reverses; i.e., reaction 1 tends to decrease and and to increase (, , ). Reaction 4 exhibits exactly the opposite behavior to reaction 1; i.e., tends to decrease and and to increase in the period and exhibits the opposite tendency in the period. Although the influence of these reactions in the period can be predicted from the relevant reduced model in Eqs. (17–20), their influence in the period cannot be predicted, since reactions 1 and 4 are not present in the appropriate reduced model in Eqs. (23–26).

Another feature that cannot be detected by inspecting the existing models relates to the influence of reactions 3f and 3b. Figure 6 shows that these reactions are shown to exercise a certain influence in the evolution of the system during the last stage of the period and throughout the one, in contrast to what is suggested by the reduced models derived in “Conventional reduced models” Section, in which these reactions are absent.

In order to validate these findings and explain the discrepancies between the reduced models in “Conventional reduced models” Section and the Importance Index findings in Fig. 6, the influence of reaction 1 and 3f in driving the system was investigated, by doubling their reaction rate constant from a specific point in time until the end of the process and by comparing the perturbed concentration profiles with the unperturbed ones. In general, a perturbation - such as the one considered here - will influence both the constraints and the manner in which the process evolves within these constraints. Naturally, the response of the constraints will manifest in accordance to the fast time scales that characterized their formation, while the response of the evolution within the constraints will be will manifest according to the slow time scales. The results such a comparison are shown in:

Figure 7, where the effect of the perturbation imposed by doubling the rate constant of , enforced from (where ) is displayed in the top row and that enforced from (where ) is displayed in the bottom row and

Figure 7

Concentration profiles for the reference case (parameter values in Table 1, solid curves) and for the case in which the reaction rate constant of is doubled (dashed curves). The perturbation is enforced from (top) and (bottom). The indicated and periods refer to the reference case.

Figure 8, where the effect of the perturbation imposed by doubling the rate constant of , enforced from (where ) is displayed.

Figure 8

Concentration profiles for the reference case (parameter values in Table 1, solid curves) and for the case in which the reaction rate constant of is doubled (dashed curves). The perturbation is enforced from , as indicated by the vertical line. The indicated and periods refer to the reference case.

Perturbation on

The top row of Fig. 7 shows that when the perturbation on is initiated at (where ), the perturbed profile of exhibits a sudden increase, while those of the other three variables exhibit a smooth departure, relative to the unperturbed profiles. The sudden increase of is due to the response to the perturbation of the equilibrium in Eq. (37) that develops in the period:where in the unperturbed case and in the perturbed one. According to the values displayed in Fig. 5, the contribution of in this equilibrium is relatively small, but not negligible; i.e., at when the perturbation is enacted , and . As a result, a perturbation on will produce a weak response. In addition, according to the values displayed in Fig. 5, the pointed variable for this equilibrium is ; . Therefore, it is this variable that will adjust the most as soon as the equilibrium in Eq. (45) is perturbed, exactly as displayed in the top row of Fig. (7); i.e., a sudden increase of produces an equally fast increase of , which is limited in size, given the contribution of to the equilibrium in Eq. (45). The displayed in Fig. 7 smooth increased rates of consumption of and and of production of and (the latter after the sudden increase due to the perturbed equilibrium), are all in agreement with the Importance Index findings reported in Fig. 6; i.e., , , , .

The bottom row of Fig. 7 shows that when the doubling of the reaction rate constant of is initiated at (), the perturbed profile of , and exhibit a sudden change and then a smooth variation, while the profile of exhibits a smooth departure, relative to the unperturbed profiles. The sudden changes are due to the the two equilibria and in Eq. (44) that develop in the period:where in the unperturbed case, in the perturbed one and denote smaller contributions. According to the values displayed in Fig. 5, the contribution of in the second equilibrium is significant; i.e., at when the perturbation is enacted , , , and . In addition, according to the values displayed in Fig. 5 and Table 4, for the second equilibrium there two the pointed variables at the start of the period, and , and only one during the remaining period, . Given that the perturbation is enacted right after the start of the period, a sudden increase of will result in a sudden increase of , via the second equilibrium in Eq. (46) and a sudden decrease of via the conservation law . In turn, this sudden increase of will result to a sudden increase of via the first equilibrium in Eq. (46), given that this is the pointed variable, as shown in Fig. 5 and Table 4. The smooth increased rates of consumption of , and and production of , exhibited either right after the perturbation (applies to ) or right after the sudden response to the perturbation of the two equilibria (applies to , and ) are in full agreement with the Importance Index findings reported in Fig. 6; i.e., , , and .

Figure 8 displays the response of the concentration profiles to the doubling of the reaction rate constant of , enacted at (where ). It is shown that the perturbation generates a sudden response of , which is followed by a slower response of , and and a subsequent even slower response of all four variables. These three successive types of responses are due to the fact that at the time when the perturbation is enacted there are two different equilibria established, which are characterized by two different time scales (), and that the perturbation modifies both of them. In particular, in the period the established equilibria and , are expressed as:so that now associates to ( in the unperturbed case, in the perturbed one) and is present in the expression of both equilibria. According to the and values displayed in Fig. 5, the contribution of in the first equilibrium in Eq. (47) is more significant than the second one in Eq. (48); i.e., at when the perturbation is enacted and . Therefore, the response to the perturbation considered here of the first equilibrium will be more pronounced than that of the second; both in the speed of its manifestation (since ) and in magnitude (since ). These features are displayed Fig. 8, where an increase of is shown to generate initially a very fast decrease of (the variable exhibiting the largest values, according to Fig. 5) via the first equilibrium and a slower decrease of and increase of (the variables exhibiting the largest values, according to Fig. 5) via the second equilibrium and the conservation law . The evolution of the four variables after these two initial transients are in accordance to the Importance Index results in Fig. 6. Specifically, the consumption of , , and the production of slow down (, , and ).

The results presented in “The M = 1 period” to “The reactions driving the process” Sections demonstrated that CSP can algorithmically provide all information needed in order to acquire system-level understanding. CSP provided everything the cumbersome analytical methodologies might provide and many more. Specifically: The significance of the identifications in item (1) is profound, especially when the mathematical model of interest is large and complex. Item (2) addresses a feature that the conventional methodologies cannot handle; that is, the fact that reactions have usually a fast and a slow component, each of which might exercise a different influence. This feature was first encountered in the case of the glycolysis cycle, where the slow component of certain reactions were shown to promote a larger cycle, while their fast component were opposing this[39]. It was shown that an increase of the reaction rate constant of these reactions was promoting their fast component, reducing significantly the period of the cycle. Similar features are manifested in the simple PK model considered here, resulting to (1) the reversal of the influence of reactions 1 and 4 during the transition from the to the period and (2) the influence in the period of reaction 3, which was deemed equilibrated there; .

The analytical methodology does not provide a method for the identification of the fast variables and fast reactions. These identifications are usually based on the experience and intuition of the investigator. In contrast, via the Partitipation Index , the Pointer and Time scale Participation Index , CSP can make these identifications algorithmically, as it was demonstrated in “The M = 1 period” and “The M = 2 period” Sections.

The analytic methodology assumes that a reaction deemed fast must be eliminated from the reduced model. As a result its influence on the slow evolution cannot be assessed. In contrast, CSP allows for a reaction to have a fast and a slow component, so its influence in shaping the constraints and in driving the system within these constraints can be assessed with the Importance Index , as it was demonstrated in “Perturbation on ” Section.

Specifically, according to Fig. 6, towards the end of the period reaction 1 tends to decrease and to increase and , while in the period this influence of reaction 1 reverses. A similar behavior is exhibited by reaction 4; i.e., in the period reaction 4 tends to increase and to decrease and , while it tends to the opposite in the period. This feature is due to the fact that reactions 1 and 4 have a fast and a slow component, which tend to move the system in opposing directions. In particular, Fig. 6 shows that in the and periods reaction 1 contributes significantly to both and , while has negligible contribution to the fastest amplitude ; as indicated by the large values of and and the negligible ones . This implies that reaction 1 has non-negligible components only in the 2nd and 3rd CSP modes:where the first term in the RHS of this expression belongs to the relatively fast 2nd CSP mode () and the second term belongs to the relatively slow 3rd one (). In the period, according to Eq. (29) the system is formally driven by both these CSP modes. However, according to the amplitude profiles in Fig. 4, , so in this period it is mode that dominates. Therefore, the component of reaction 1 in the relatively fast 2nd mode exhibits much largest influence in the period, in comparison to its component in the relaatively slow 3rd mode. As a result, the value of the Importance Index in this period reflects the influence of the fast component of reaction 1. In the period the 2nd CSP mode is exhausted (i.e., ) and the system is driven only by the 3rd mode , in accordance to Eq. (30). As a result, the influence of reaction 1 is exercised only through its component in the relatively slow 3rd mode and there reflects the influence of the slow component of reaction 1. The fact that the values of , and have opposite signs in the and periods, as displayed in Fig. 6, demonstrate that the fast and slow components of reaction 1 exhibit opposing influences. A similar explanation can be provided for the opposing influence of reaction 4 in the and periods, as shown in Fig. 6.

Reaction 3 exercises a significant influence in all three modes during the period, as it is demonstrated by the non-negligible Amplitude Participation Indices and that are displayed in Fig. 5 for . This indicates that this reaction has non-negligible components in all three modes:where . The values of and in Fig. 5 for the fastest mode during the period indicate that this reaction is in partial equilibrium; . However, this does not preclude this reaction having an influence in the slower modes, as it it demonstrated by the non-negligible Importance Indices and , as shown in Fig. 6 and in more details in Fig. 8.

The traditional methods cannot distinguish the fast and slow components of the reactions. Therefore, the possibility that a change of the value of a parameter to have a certain influence in the short-run and the opposite in the long-run cannot be discovered. In contrast, the algorithmic CSP methodology can identify both (1) the component of a reaction that associates to a certain time scale via the Partitipation Index , the Pointer and Time scale Participation Index and (2) the influence of this component in driving the process via the Amplitude Participation Indices .

Dominant dynamics in the M = 1 and M = 2 periods

It was demonstrated in the previous sections that CSP was able to algorithmically (1) identify the constraints that develop due to the fast dynamics and (2) recognize their nature; i.e., the fast variable, the fast reaction and the reactions participating in the constraints. It was also demonstrated that CSP can provide very accurate expressions for the constraints, which is reflected in the accuracy of the reduced model; see Eq. (31) for the general form of the constraints and the reduced model.

It will be demonstrated next that the CSP generated diagnostics, displayed in Figs. 2, 3, 4 and 5 and Tables 2 and 4, can provide the required knowledge in order to predict the response of the system when subjected to a change in the value of a parameter. The demonstration will address first the period and then the one.

Let us first consider the period, for which is was shown that the process evolves within one constraint in the form of Eq. (35):For this case, the Amplitude Participation Index displayed in Figure 5 shows that only reactions 1, 3f and 3b provide significant contributions in the occurring cancellations, i.e.:where , since and , and throughout the period; see Table 2. As the Pointer diagnostics reported in Figure 5 show, the pointed variable for this equilibrium is . Therefore it is this variable that will adjust the most in this period when a value of a parameter of those appearing in Eq. (52) is perturbed. For example, an increase of () will cause an increase of , while an increase of () will have no effect, since it is not present in this equilibrium, Eq. (52). These conclusions are validated by the results displayed in Figs. 9 and 10, where concentration profiles computed with the reference parameter values in Table 1 are compared with profiles computed with perturbed values and ). It is shown that variations of these two parameters have no other significant effect in the period where , other than the increase of in the case of the perturbed value .

Figure 9

The reference (solid, ) and the perturbed (dashed, ) concentration and time scale profiles. The indicated and periods refer to the reference case.

Figure 10

The reference (solid, ) and the perturbed (dashed, ) concentration and time scale profiles. The indicated and periods refer to the reference case.

Regarding the influence of the perturbations to the time scales (1) that characterizes the development of the 1st equilibrium and (2) that characterizes the evolution of the system in the period, the changes in their profiles displayed in Figs. 9 and 10 are in full agreement with:

The Time scale Participation Index diagnostics for reported in Fig. 5, according to which the related to the equilibrium in Eq. (52) time scale depends only on reaction 3f; i.e., is independent of or . Therefore, there will be no change in the period on how fast the first equilibrium will be established when or are perturbed, as shown in Figs. 9 and 10.

The Time scale Participation Index diagnostics for reported in Fig. 5, according to which the time scale , which characterizes the second mode that drives the system in the period (since ), depends mainly on reaction 4; i.e., is also independent of or . Therefore, there will be no change on how fast the system evolves in the period or on the time required for the second equilibrium to be established when or are perturbed, as shown in Figs. 9 and 10.

Let us next consider the period, for which is was shown that the process evolves within two constraints in the form of Eqs. (35) and (41):The Amplitude Participation Index displayed in Fig. 5 shows that in the period considered here only reactions 3f and 3b provide significant contributions in the occurring cancellations in the first equilibrium (, Eq. (44)) and only reactions 1, 2 and 4 contribute to the second (, Eq. (44)), i.e.:where and throughout the period and ; see Table 4. As the Pointers and reported in Fig. 5 indicate, the pointed variable for the 1st equilibrium throughout the period is and for the 2nd equilibrium is mainly and at the start of the period as well. Therefore, these are the variables that will adjust the most in this period when a value of a parameter is perturbed.

Considering the increase of (), Eq. (55) shows that this will cause an increase of the pointed variable via the 1st equilibrium. This increase will cause a decrease of via the 2nd equilibrium Eq. (56), which in turn will cause an increase of via the conservation law . These CSP-based predictions are in full agreement with the findings displayed in Figure 9.

Considering now the increase of (), Eq. (56) shows that this will cause a decrease of the pointed variable via the 2nd equilibrium, which in turn will cause an increase of via the conservation law . The decrease of and the increase of will in turn cause a decrease of the pointed variable via the 1st equilibrium Eq. (55). Again, these CSP-based predictions are in full agreement with the findings displayed in Figure 10.

Note that in both perturbations considered, the profile follows that of . This feature is in accordance to the 2nd equilibrium Eq. (56), in which the contribution of reaction 2 decreases as the process moves deeper into the period; i.e., , as shown in Fig. 5.

Regarding the influence of the perturbations to the time scales (1) and that characterize the development of the 1st and 2nd equilibria and (2) that characterizes the evolution of the system in the period, the changes in their profiles displayed in Figs. 9 and 10 are in full agreement with: These predictions are validated by the results shown in Figs. 9 and 10. In particular, it is shown in these figures the fastest consumption of , and and the fastest reach to steady state of in the period, which are due to the faster characteristic time scale .

The Time scale Participation Indices and diagnostics reported in Figure 5, according to which the related to the first equilibrium time scale depends only on reaction 3f and the related to the second equilibrium time scale depends basically on reaction 4 (i.e., both are independent of or ).

The Time scale Participation Index diagnostics reported in Figure 5, according to which the time scale that characterizes the only active mode depends mainly on all reactions; i.e., it depends on and . Since and are both negative, increasing or (, ) will result in a faster ; i.e., the process will evolve faster in the period.

Conclusions

A demonstration is provided here on how algotithmic multi-scale analysis can provide system-level understanding. This work refers to the traditional analytic and to an algorithmic (CSP) treatment of a simple model that exhibits multi-scale dynamics. It is shown that CSP provides algorithmically all the results that are provided by analytical techniques along with additional ones, that cannot be obtained analytically. In contrast to the traditional analytical treatment, CSP was shown to require no input or experience from the investigator.

In a multi-scale process, the fast and slow dynamics exhibit distinct influence: (1) the fast dynamics generates constraints in which the process evolves and (2) this evolution is governed by a slow model (free of the fast dynamics). A perturbation along the fast components of the model (reactions) will have an immediate effect that will affect the fast variables, since it will modify the existing constraints, the development of which is characterized by the fast time scales. In contrast, a perturbation along the slow components will have a gradual effect and can potentially affect all variables, since it will modify the evolution of the process within the established constraints. A multi-scale analysis aims in identifying (1) the constraints, (2) the slow system that governs the evolution within the constraints, (3) the reactions whose fast component is responsible for the development of the constraints, (4) the reactions whose slow component determines the slow evolution, (5) the reactions that can be ignored; i.e., those that do not contribute to the fast or slow dynamics and (6) the fast variables; i.e., those that are “slaved” to the rest. CSP provides all these in an algorithmic fashion, so no input from the investigator is required and the analysis can be repeated for any parameter set of interest. The conclusions reached for a given set of parameters are valid for in a large domain of the parameter space, in which basic dynamics properties (such as time scale separation) are not altered.

It was demonstrated that the CSP diagnostics provide all necessary information to predict the response of the system when a change in the parameters is introduced with the use of specific algorithmic tools; (identifies the reactions participating in the n-th equilibrium), (identifies the variables that will adjust the most when the system is perturbed), (identifies the reactions that generate the n-th time scale) and (identifies the reactions that drive the n-th variable within the established equilibria). These numerical diagnostics can provide all required information in order to acquire system-level understanding and seek ways to redesign the process.

Finally, it was shown that CSP can distinguish the dual role a reaction might exhibit, via its opposing influence in the fast and slow dynamics. This is a feature that popular methodologies (such as, paper-and-pencil singular perturbation analysis, QSSA and PEA) cannot handle, since they cannot consider reactions having a significant fast and slow component. The same applies to sensitivity analysis, which cannot distinguish fast and slow reactions.

It is important to note that our automatic algorithmic analysis came to some conclusions that are very similar to the ones a pharmacokineticist would have derived based on detailed manual analysis of the system. For example, in the very early time points after antibody administration (i.e., period) the entry of antibody from plasma to endosome via pinocytosis (i.e., reaction 1) affects plasma PK the most, and minimizing this process with help achieve higher plasma concentrations. This aspect can be leveraged to design better antibodies for more efficient elimination of pathogenic targets. For example, recently sweeping antibodies have been developed that utilize surface FcRn or Fc-gamma receptors to increase cellular uptake of antibodies and associated pathogenic protein molecules (e.g., PCSK9, C5 etc.)[61,62]. Literature and our lab’s recent data, which show that positively charged antibody molecules demonstrate a higher cellular uptake rate and a reduced plasma exposure, also support the conclusion derived by the model over here[63,64]. As the time progresses (i.e., period), reducing the pinocytosis in the cells and increasing the recycling of FcRn bound mAbs (i.e., strengthening reaction 4) both will help maintain higher plasma exposure of mAbs. Even at later time points (i.e., period), weakening reaction 1, strengthening reaction 4, weakening reaction 2 (i.e., reducing lysosomal degradation rate of mAbs), and strengthening reaction 3f along with weakening reaction 3b, will all help improve the plasma exposure of mAbs. This aspect of the process has been heavily utilized by the protein engineers to develop antibodies with very long half-lives. For example, mutations (e.g., YTE, LS etc.)[46,65] have been introduced in the Fc region of the antibodies to improve their binding to the FcRn at pH 6.0, leading to half-life of several months (compared to 3 weeks for wild-type antibodies in human). Our analysis also provides additional insight not generated by other methods. For example, identification of the period in which the influence of each reaction in increasing plasma PK of mAb is exercised (e.g., reaction 4 becomes influential when enough FcRn bound mAb complex has been accumulated). These unprecedented insights can further stimulate protein engineering efforts to discover and develop better antibody-based drug molecules.

This work demonstrates the ability of CSP in fully analyzing complex Pharmacokinetic and Quantitative Systems Pharmacology models. Due to the algorithmic nature of CSP, the type of analysis introduced here can be extended to much more complex models; i.e., models that include target mediated drug disposition (TMDD) to characterize the PK of target binding antibodies and antibody drug conjugates (ADC)[66-69]. The new features introduced by these models will be represented as additional terms in the governing Eq. (1). The influence of these terms can be easily assessed by the algorithmic tools of CSP.

Supplementary Information.