A mathematical analysis of rebound in a target-mediated drug disposition model: II. With feedback.

We consider the possibility of free receptor (antigen/cytokine) levels rebounding to higher than the baseline level after the application of an antibody drug using a target-mediated drug disposition model. It is assumed that the receptor synthesis rate experiences homeostatic feedback from the receptor levels. It is shown for a very fast feedback response, that the occurrence of rebound is determined by the ratio of the elimination rates, in a very similar way as for no feedback. However, for a slow feedback response, there will always be rebound. This result is illustrated with an example involving the drug efalizumab for patients with psoriasis. It is shown that slow feedback can be a plausible explanation for the observed rebound in this example.

Introduction

In this paper, we continue our investigation into the phenomenon of receptor rebound, i.e., a post-dose rise in receptor levels to higher than pre-dose (baseline). In our previous paper (Aston et al. 2014), we showed that if no homeostatic feedback is present, rebound will occur if and only if the elimination rate of the target-drug complex is slower than both the elimination rate of the drug and the elimination rate of the target. Binding to cell-surface receptors typically results in accelerated turnover of the anti-body and hence complex, so it can be expected that rebound will occur in rare cases only. However,  Ng et al. (2005) describe a treatment for psoriasis patients with efalizumab in which rebound is observed. They also derive a model that shows good agreement with the experimental observations. A main difference between this model and standard TMDD models is that the model includes receptor feedback for the synthesis rate. In the basic TMDD model, reduction of free target levels does not have any impact on the endogenous production or elimination rate of the free target, i.e., no endogenous feedback control exists to compensate for the antibody effect on target. The model in Ng et al. (2005) includes such feedback control via an additional differential equation for the synthesis rate. Another indicator of the importance of feedback in the synthesis rate is the recent paper by Kristensen et al. (2013), in which it is postulated that the protein synthesis rate is the predominant regulator of protein expression during differentiation.

In this paper we investigate how receptor feedback influences the occurrence of rebound in the TMDD model. The receptor feedback is usually a dynamical process via some moderators and we modify the TMDD model to include feedback by adding an additional differential equation. If the feedback is very fast, a quasi-equilibrium approach can be used and the feedback can be included in the synthesis term itself Hek (2010). This approach is similar to the one leading to the Michaelis–Menten approximation or the Quasi-Steady-State approximation. We will use the term “direct feedback” for the quasi-equilibrium approximation. The modification to include feedback is presented in full detail in Sect. 3, after a review of the main results of our first paper (Aston et al. 2014) on rebound in the basic TMDD model without feedback in Sect. 2. A wide class of feedback functions is considered with some natural assumptions on their action.

The TMDD model with the direct feedback approximation is analysed in Sect. 4. This analysis shows that the existence or non-existence of rebound is still linked to the elimination rates, though rebound can be expected in a larger region in the elimination parameter plane, see the left plot in Fig. 1. The general case of feedback via a moderator is analysed in Sect. 5. Now the response speed of the feedback moderator plays an important role as well as the elimination rates. It is shown that there is a similar region in the elimination rate plane as for the basic TMDD model for which rebound will occur for any response speed. Furthermore, if the feedback responds slowly to a change in the receptor levels, rebound will occur for any value of the elimination rates. A schematic overview of this result is in the right plot of Fig. 1.

Fig. 1

A schematic overview of the elimination plane. On the horizontal axis is the elimination rate of the ligand (antibody/drug), denoted by , and on the vertical axis is the elimination rate of the receptor (antigen/target), denoted by . The elimination rate of the antibody–antigen complex (drug–target product) is denoted by . On the left is an overview of the occurrence of rebound in the case of the direct feedback approximation. In the red region, rebound will occur and in the green region no rebound can occur. The symbols h, m and M are related to the feedback function and will be defined in Sect. 4. In the case where there is no feedback, and the region in which there is rebound is the red region above the line . On the right is an overview of the occurrence of feedback, depending on the speed of the feedback response. There is a red region in the elimination plane, that equals the direct feedback rebound region, for which rebound will occur for any speed of the feedback response. In addition, if the response is slow, rebound will occur for any elimination values. The precise response speed for which the rebound stops occurring will vary across the pink region (color figure online)

The TMDD model is a one-compartment model. In Sect. 6, we consider briefly two classes of more general models with feedback (one class is related to multiple compartments, the other to more general feedback mechanisms) and show that rebound will generically occur in such models for slow feedback moderator response. In this paper, the word generic refers to the assumptions that the linearisation about the baseline state has no degenerate eigenvalues and that trajectories approaching the attracting baseline state do so tangent to the eigenvector associated with the least attracting eigenvalue. Section 7 illustrates the theory for two examples. The first example is a standard TMDD model for the IgE mAb omalizumab. We discuss what the effect of the two types of feedback mechanisms would be in this case . The second example is the model in Ng et al. (2005) which describes the efficacy of efalizumab for patients with psoriasis. Simulations show that the model does not lead to rebound if the feedback is turned off, but a significant rebound (about 140% of baseline) will occur if feedback is included. These results can be predicted with our analysis and underpin the observations of rebound in some patients. Finally, Sect. 8 contains a discussion of the results obtained in this paper and poses some open questions. The proofs of some of the technical results in the paper are given in the Appendix.

Review of previous results

In our previous paper (Aston et al. 2014), we considered a one-compartment TMDD model based on the original work of Levy (1994) where the ligand L binds reversibly with the receptor R to form a receptor-ligand complex P as shown in Fig. 2. The TMDD model assumes a mechanism-based reaction to explain the ligand-receptor interaction. The parameters of the model are the binding rate constants and , the receptor turnover and elimination rates and , and the elimination rates of the ligand and complex and . The system is assumed to be initially at steady state, into which a single bolus infusion of the ligand into the central (plasma) compartment is made (represented in Fig. 2 by ‘In’). The differential equations that comprise the mathematical model for this system are given by A steady state of this system is given by , . Adding the bolus injection gives the initial conditionsAfter the ligand is added to the system in its baseline state, initially the receptor level decreases, but after a while it goes up again and returns to its baseline value, since the steady state is globally asymptotically stable (Aston et al. 2014). Rebound occurs if, in the return to the baseline value, the receptor level increases to values above the baseline . In Aston et al. (2014), we determined precise conditions for the existence and non-existence of rebound. Our main result was the following.

Fig. 2

The TMDD reaction mechanism

Theorem 2.1

Rebound occurs in Eqs. (1)–(3) if and only if

This result shows that rebound occurs if and only if the elimination rate of the product is slower than the elimination rates of both the ligand and the receptor. This is represented graphically by different regions in the plane, as shown in Fig. 3. The above condition for rebound in terms of non-dimensional parameters is given bywhere

Fig. 3

The green region shows the part of the parameter plane where there is no rebound, and the red region shows where there is rebound. The boundaries are given by the lines and . This figure is similar to Fig. 2 in Aston et al. (2014) (color figure online)

Intuitively, if the clearance of the complex is slow relative to the other two clearance rates, when there is a significant build-up of the complex due to the rapid binding of ligand and receptor, then this complex is only eliminated slowly. This results in an additional “production” of the receptor when the complex dissociates which is likely to lead to rebound. Our analysis quantifies this and we now show how these thresholds that we derived previously are affected if receptor feedback is added to the system.

Including feedback in the TMDD model

We now extend our previous work on rebound in the basic TMDD model by extending the model to include feedback on the synthesis rate. The feedback gives an increase in the production rate of the receptor if the receptor levels drop under the baseline. Similarly, if the receptor level exceeds the baseline value, then the feedback induces a decrease in the production rate. Thus, this mechanism can be classified as negative feedback.

We start with the TMDD equations (1)–(3) but now assume that the rate of production of R () is not constant but varies as R changes. Thus, Eq. (2) is now modified to have the formThe new feedback variable F, which will be related to R, has the effect of varying the production rate of the receptor. The feedback is usually via some moderator and its dynamics is described with the differential equationThe constant determines the strength of the feedback and determines the time scale on which the feedback moderator responds to changes in the receptor levels. Clearly, if , then F is a constant and the model reduces to the basic TMDD model with no feedback. If is small, the feedback response is slow, while if is larger, there is a very fast feedback response. Formally, dividing (6) through by and taking the limit , we see that the derivative term vanishes and the fixed relationis obtained. This motivates a quasi-equilibrium approximation in which the differential equation for F is replaced with the relation (7). We will call this the direct feedback approximation.

Next we consider the function H(R). The function must be such that the model has the following features:This leads to the following general assumptions on the function H(R):We note that assumption (H1) ensures that in case of the direct feedback approximation (i.e, F is a function of R and (7) holds) the feedback , so that the production rate of R is always positive. Alternatively, if the dynamics of F is described by the differential equation (6), then this condition ensures thatand so the line cannot be crossed from above to below. This also ensures that for all and hence the production rate is always positive.

the production rate of the receptor, denoted by , should be positive;

the receptor concentration should have the baseline ;

the feedback mechanism should reduce the production rate of the receptor when the receptor level exceeds the baseline and should increase it when the receptor level is below the baseline.

for all ;

, where ;

when and when ;

.

The basic TMMD model without feedback, Eqs. (1)–(3), has a steady state solution given by , . This will still be a steady state of our modified equations in both the full model and the direct feedback approximation, provided that (H2) holds. For the full model, we observe that the steady state value for F is . Finally, the assumption (H3) ensures that the production of the receptor slows down when there is excess receptor, and increases when the receptor level is below the baseline and (H4) is a technical assumption assuring sufficient smoothness. Combining (H2)–(H4) implies that . We note that if H is a positive, strictly decreasing, smooth function with , then (H1)–(H4) are satisfied. However, the stated conditions are more general and do not require H to be strictly decreasing.

We now have two modified TMDD models that incorporate feedback. For the direct feedback approximation, we use (7) to replace F with H(R) which gives the differential equations Alternatively, the full feedback model involves an extra differential equation for the moderator F and the modified TMDD model is given by The initial conditions for both models are given by (4), with the extra condition for the feedback moderator.

We non-dimensionalise these equations in the same way as in Aston et al. (2014). We note that our new variable F is non-dimensional, and so we define the dimensionless variablesWe also define a new functionThe above assumptions on the function H(R) translate into the following assumptions on h(y):For later use, we defineand note that the assumptions (h2)–(h4) imply that .

for all ;

;

when and when ;

.

In terms of the non-dimensional quantities, the direct feedback equations (8)–(10) become with initial conditionswhere dot denotes differentiation with respect to and the dimensionless parameters are defined as in Aston et al. (2014) byClearly, this choice of non-dimensionalisation requires that and (so that ). We note that all parameters must be non-negative due to their physical meaning, and additionally we will assume that they are in fact all strictly positive. The three variables x, y, and z are related to physical quantities and so must also be non-negative.

Similarly, the non-dimensional equations for the full feedback model (11)–(14) are given by with initial conditionsFurthermore,is a non-dimensional measure for the response speed of the feedback moderator.

Analysis of the TMDD model with direct feedback

In this section, we will analyse the TMDD model with the direct feedback approximation, i.e., the Eqs. (16)–(18), with initial conditions (19), and seek to determine conditions for the existence or non-existence of rebound in this model. However, we start by proving some basic results for the model equations.

Invariance, steady states and stability

For the basic TMDD model, we proved in Aston et al. (2014) that the positive octant is invariant; that there is a unique steady state in the positive octant; and that this steady state is a global attractor. In this section we show that the same properties hold for the TMDD model with direct feedback. The proofs of these results can be found in Appendix A. We start by stating invariance of the positive octant.

Lemma 4.1

The octant of defined by is invariant under the flow of Eqs. (16)–(18).

This result shows that the equations are a good model in the sense that the concentrations of the ligand, complex and receptor can never go negative. Next we consider the steady states of Eqs. (16)–(18).

Lemma 4.2

In the region of defined by , the Eqs. (16)–(18) have a unique steady state given byThis steady state is globally asymptotically stable.

This result shows that for all non-negative initial conditions and for all positive parameter values, each of the variables will converge to their unique steady state value.

It is also straightforward to show that the eigenvalues of the Jacobian of Eqs. (16)–(18) evaluated at the steady state (26) are very similar to those given in Aston et al. (2014). In particular, the two eigenvalues and are the same as previously and are given by while the third eigenvalue, which previously was , is now given bywhere we recall that . The corresponding eigenvectors are given by

Receptor dynamics

Before considering the full model, we first consider the effect of feedback on the dynamics of y, the non-dimensional form of the receptor, in the absence of the ligand or product. As with the basic TMDD model, we note that satisfies (16) and (17). Substituting into (18), we obtain the single differential equationAssumption (h2) ensures that the previous steady state value is also a steady state of this equation while assumption (h3) ensures that this is the only non-negative steady state solution and that if and if . The theory of scalar, autonomous differential equations (Jordan and Smith 2007) then tells us that if the initial condition is less than the steady state () then will increase monotonically to the steady state, whereas if is greater than the steady state (), then will decrease monotonically to the steady state. As is monotonic, there is no possibility of rebound occurring in the dynamics of the receptor on its own.

Assumptions (h2) and (h3) together with definition (15) imply thatAs , the inclusion of the function h(y) ensures that the asymptotic rate of convergence to the steady state (proportional to ) will be faster than or equal to the rate of convergence which occurs when (proportional to ). Thus, qualitatively, the direct feedback has the effect of speeding up the monotonic convergence to the steady state in the absence of ligand and product, but there is no rebound in this case, i.e., if , then for all .

Conditions for rebound

Having established that feedback does not lead to rebound in the absence of the ligand, we now investigate the effect of adding the ligand. A local approximation of the feedback near the baseline will give a linear feedback function. So we first consider the special case that the feedback function h(y) is linear for . It turns out to be straightforward to extend our results for the TMDD model without feedback to this case. The general nonlinear case will be analysed next and builds on the ideas of the linear analysis.

A mainly linear feedback function

We consider a feedback function which is linear on a y-interval including [0, 1]:where , , and , , withThis function satisfies all the conditions (h1)–(h4). For the linear part of the function h, we note that and so Eq. (18) becomesWe then have the following result.

Theorem 4.3

Rebound occurs in Eqs. (16), (17), (33) if and only if

Proof

We note that for , Eqs. (16), (17) and (33) are the same as the basic TMDD model that we studied in Aston et al. (2014) except that has been replaced by . We claim that the results of Theorem 3.1 of Aston et al. (2014) apply in this case also after making this parameter replacement. To see this, we note that when proving the absence of rebound the only relevant values of y are and so the matching point in (33) is not reached and the feedback function for is therefore not relevant. Moreover, the proof of rebound in Aston et al. (2014) involves only the asymptotic behaviour near to the globally stable steady state and the dynamics for will give all the necessary information to determine this asymptotic behaviour. Hence, we can conclude that Theorem 3.1 in Aston et al. (2014) can be applied in this situation with replaced by .

Converting back to the dimensional parameters gives the following result.

Corollary 4.4

Consider the mainly linear feedback functionwhere , , , and and its first and second derivatives match the linear function at .

Rebound occurs in the model given by Eqs. (8)–(10) with this feedback function if and only if

The rebound region in this case is shown in Fig. 4. We note by comparison with Fig. 3 that this mainly linear feedback function has enlarged the region in which rebound occurs compared to the case with no feedback.

Fig. 4

The green region shows the part of the parameter plane where there is no rebound, and the red region shows where there is rebound for H(R) given by (34). The solid line denotes the boundary of the rebound region when no feedback is present and the dashed line denotes the boundary of the rebound region when feedback is present (color figure online)

A general nonlinear feedback function

We now consider the more general case where the feedback function h is nonlinear. We express the function h(y) in the formDifferentiating (35), evaluating at and using (15), we see thatCondition (h3) implies thatWith this form of h, Eq. (18) becomes

Remark

We note that expanding h(y) in a Taylor series with remainder gives that where lies between y and 1. However, the fact that is an unknown function is not helpful, which is why we use the form of h given in (35).

We now consider regions similar to the ones defined for the model without feedback in Aston et al. (2014). First we consider region II and show that no rebound occurs in this region.

Theorem 4.5

Rebound does not occur if .

Following through the proof of Theorem 3.4 of Aston et al. (2014) for these modified equations, the only difference is thatNow by (37) and so when , which is the condition required on this plane in the proof of Theorem 3.4 of Aston et al. (2014). Thus, the result holds in this case also.

The next result considers a region similar to region III and we show that rebound does occur.

Theorem 4.6

Rebound occurs if and , where is given by (27).

The Jacobian matrix for Eqs. (16), (17) and (38) evaluated at the steady state is the same as that for the simple TMDD model given in Aston et al. (2014) but with replaced by . The proof of Theorem 3.5 of Aston et al. (2014) (region III) is entirely based on the eigenvalues and eigenvectors of this matrix, and so the result holds in this case also, with the appropriate replacement of . Since does not involve , there is no change to this eigenvalue.

The proofs in regions I and IV of the parameter plane in Aston et al. (2014) both made use of the fact that for the model without feedback the plane is invariant when . However, this is no longer the case for our modified equations which include feedback. The differential equation for the dimensionless variable associated with the total amount of receptor is now given byand soClearly no condition on the parameters will give this derivative to be zero (unless is constant, the case considered in the previous section), and so we see, as already stated, that the plane is no longer invariant. However, the function h(y) has a continuous derivative by (h4) and this implies that is continuous. Thus the Extreme Value Theorem for continuous functions gives that has a minimum and a maximum for . So we defineUsing (36), this immediately implies thatAlso, since for all by (37), we getSince m and M are the minimum and maximum of , then clearlyand substituting for from (35) and rearranging givesGraphically, this means that on the interval , the function h(y) is bounded by two straight lines, both of which pass through the point (1, 1), and which have slopes and . This situation is illustrated in Fig. 5. Note that the Mean Value Theorem gives that and .

Fig. 5

The function h(y) (solid curve) is bounded by two straight lines (dashed lines) with slopes and for

With the bounds m and M, we can give a slightly weaker version of Lemma 2.3 in Aston et al. (2014) and get regions in the -plane for which the dimensionless variable v associated with the total amount of receptor stays on one side of the plane .

Lemma 4.7

For all , the dimensionless variable v associated with the total amount of receptor satisfies

We consider the vector field at the plane . Equation (40) givesSince we have restricted attention to the plane , this implies that and so as z must be non-negative. As is non-negative, we can use the bounds on given in (41) to getFirst we consider the case . In this case, (45) gives , hence the vector field on the plane is either tangent to or points towards . This implies that the region of phase space defined by , is invariant in forward time since the vector field on the four boundary planes is either tangent to the plane or points into the region (Smith 1995). This was proved for the planes , and in Lemma 4.1 and for the plane above. Since and this region is invariant, we conclude that for all .

Next we consider the case . Now the other inequality in (45) gives , hence the vector field on the plane is either tangent to or points towards . This implies that the region of phase space defined by , is invariant in forward time. Since , this implies that for all .

With this observation, we can now derive a result on rebound for regions similar to regions I and IV in Aston et al. (2014).

Theorem 4.8

Define m and M as in (41), i.e., the function h(y) satisfies (44).

Rebound does not occur if

Rebound does occur if

Using (42) and (43), we note thatThis illustrates in particular that there is no overlap between the region for which we have shown the existence of rebound and the region for which we have shown that there is no rebound.

Proof of Theorem 4.8

The case when (46) holds is similar to region I in Aston et al. (2014) and has the property that for all by Lemma 4.7. Thus for all , we have and so rebound cannot occur in this case.

When (47) holds (a condition similar to the one for region IV in Aston et al. (2014)), Lemma 4.7 ensures that for all . The right hand inequality in (47) implies that and so the eigenvalue closest to zero is , as occurred in region IV previously. Thus, trajectories will generically approach the steady state tangent to the y-axis, since the eigenvector associated with the eigenvalue points along this axis. As the trajectory always satisfies for all then we again conclude that generically, almost all orbits will approach from above and hence rebound occurs.

We must again eliminate the possibility that an orbit might approach the steady state tangent to one of the other eigenvectors for a particular choice of parameters. As previously, the linearised manifold in the direction of the eigenvector is outside of the phase space and so an orbit cannot approach tangent to this one-dimensional manifold.

The one-dimensional linearised manifold in the direction of the eigenvector given in (30) is We note that from the right hand inequality in (47). Thus, we must take to ensure that this linearised manifold gives positive values of x and z (since by Lemma 2.7 of Aston et al. 2014). The linearised manifold for is then given by Using (42) and the left hand inequality of (47) gives . Thus, this one-dimensional linearised manifold has a negative slope and so occurs for . Since the trajectory in this case must satisfy , it is clearly not possible for it to approach the steady state tangent to this manifold. Thus, we conclude that all trajectories must approach the steady state tangent to the y-axis from above and hence rebound will occur.

In Theorems 4.6 and 4.8, we have proved that rebound occurs on both sides of the line with . We now show that there is also rebound along this line.

Theorem 4.9

Rebound occurs if and .

The proof uses exactly the same arguments as the proof of Theorem 3.7 in Aston et al. (2014). Along the curve , the eigenvalues and collide and have one eigenvector, which is along the y-axis, and one generalised eigenvector. Thus all solutions will asymptotically align with the y-axis. Since , Lemma 4.7 gives that and the intersection of the y-axis and the region is the part of the y-axis with . Hence all orbits must approach the steady state from above and rebound will occur.

We summarise the results of Theorems 4.5, 4.6, 4.8 and 4.9 as follows.

Theorem 4.10

DefineFor the model equations (16)–(18) and for any function h(y) satisfying the assumptions (h1)–(h4)We can also express these results in the terms of the dimensional parameters. Let H(R) be a function satisfying the assumptions (H1)–(H4) and defineFor the model equations (8)–(10) with

there is no rebound if or and ;

rebound does occur if and .

there is no rebound if or and ;

rebound does occur if and .

These results are illustrated in Fig. 6.

Fig. 6

Regions in the plane where rebound does and does not occur for the nonlinear feedback function H(R). Red indicates regions where rebound will occur and green indicates regions where rebound does not occur. In the white region, rebound may or may not occur depending on the particular feedback function H (color figure online)

We now focus on a few aspects of the theorem and compare it with the case where there is no feedback (, and so ).In summary, the region where rebound occurs increases slightly to include some more values with when direct feedback is introduced into the basic TMDD model.

In the equations without feedback, rebound occurred in the region and (regions III and IV) and with feedback, rebound also occurs in this region.

In the equations without feedback, there was no rebound in the region (region II) and with feedback, there is also no rebound in this region.

In the equations without feedback, there was no rebound in the region (the part of region I that lies outside region II), but with feedback, there is rebound in some parts of this region but not in others, the precise details depend on the function h.

In the region , , we have no general results, and the existence or otherwise of rebound will depend on the particular function h.

If H is mainly linear, i.e. , for , then . In this case, the region of uncertainly disappears and we get that rebound will occur if and only if and .

If H is nonlinear but with (in which case the lower bounding curve for H(R) is the tangent to the curve at ), then the red region where there is rebound for in Fig. 4 is preserved and the region of uncertainty is below this line. Conversely, if (so that the upper bounding curve for H(R) is the tangent at ), then the green region where there is no rebound for in Fig. 4 is preserved and the region of uncertainty is above this line.

We note that when (the left green region in Fig. 6), then there is no rebound in the direct feedback model for any feedback function H(R). As discussed above, when , then the form of the feedback function determines whether or not rebound occurs.

Analysis of the full TMDD with feedback model

In the full TMMD with feedback model, we have an extra differential equation that governs the dynamics of the moderator w and so our equations are now (20)–(23) with initial condition (24). As mentioned previously, the parameter influences the time scale for the response dynamics. If is small, the feedback responds slowly and we will show that this always leads to rebound. On the other hand, if is very large, the feedback responds very fast and we will show that rebound will occur in the same region as found with the direct feedback approximation.

Dynamics of the receptor with feedback, without ligand

Before considering the full model with feedback, we first consider the effect of feedback on the dynamics of the receptor y (or R in dimensional form) in the absence of the ligand or product. We substitute into (20)–(23) to give the two equations Again, it is easily shown, using (h2) and (h3), that these equations have the unique steady state , i.e. the baseline values. Below we will prove the following rebound result.

Theorem 5.1

For any , there are always some initial conditions with which give rise to trajectories with rebound.

First we consider the linear stability of the baseline state as this gives the local behaviour near to the baseline. The Jacobian matrix evaluated at this steady state is given by(recalling that ). The eigenvalues of can be found as solutions of the quadratic characteristic polynomial derived from this matrix. The discriminant of this quadratic equation isand the eigenvalues are while the eigenvectors are of the formNote that in this section, the definition of the eigenvalue is different from that used in the previous section. The discriminant D is itself a quadratic function of which has two positive roots given byClearly the discriminant is negative between these roots and positive elsewhere, and so we have the following result.

Lemma 5.2

The Eqs. (49)–(50) have a unique steady state at and the eigenvalues and of the Jacobian matrix areFurthermore, when , and .

both real and negative if or if ;

a complex conjugate pair with negative real part if .

For a proof of the last statement, see Appendix A. Lemma 5.2 is summarised in Fig. 7 with the direction of movement of the eigenvalues indicated corresponding to increasing .

Fig. 7

Movement of the eigenvalues and in the complex plane for increasing (indicated by the direction of the arrows)

When , we see that and . As increases from zero, the two eigenvalues move towards each other and collide on the negative real axis when . After this collision, the eigenvalues become complex with , which decreases monotonically with . As increases further, the two complex eigenvalues collide on the real axis when and become real again. As increases from to , increases from the point of collision to and decreases to from this point. At the collision values , the eigenvalues are given byNote that this implies that , as .

Lemma 5.2 implies linear stability, but a stronger stability result can be proved. The proof of this result uses a Lyapunov function and is given in Appendix A.

Lemma 5.3

In the region of defined by , the steady state of Eqs. (49), (50) is globally asymptotically stable when .

Thus the trajectory for all non-negative initial conditions will eventually approach the steady state and Theorem 5.1 can be proved by analysing the solutions close to the steady states. The following lemma gives the detailed results, with the proof again in Appendix A.

Lemma 5.4

If and , then rebound in Eqs. (49), (50) will occur for all and for some values of .

If then rebound occurs for any initial condition with .

If and , then rebound in Eqs. (49), (50) will occur for some initial conditions which satisfy .

Having established the effect of the feedback on the behaviour of the receptor without ligand and product, we next consider the effect of adding the ligand.

Before considering rebound, we first state invariance and stability results for the full TMDD system with feedback, similar to the results in Sect. 4.1 for the TMDD system with direct feedback. The proofs of the statements in this section can be found in Appendix A.

The first three of our variables are related to physical quantities and so must be non-negative and we have also assumed that the feedback variable w is non-negative. Therefore, we now establish the invariance of the region with for Eqs. (20)–(23).

Lemma 5.5

The region of defined by is invariant under the flow of Eqs. (20)–(23) when .

As with our previous model, this again shows that none of the variables can go negative, which is an essential property of a good model. We now consider the steady states of Eqs. (20)–(23).

Lemma 5.6

In the region of defined by , there is a unique steady state of Eqs. (20)–(23), given by

The Jacobian matrix obtained from Eqs. (20)–(23) evaluated at the steady state solution (57) is given by(recalling that ). The eigenvalues of this matrix can be found from the top left matrix and the bottom right matrix. From the top left matrix, we obtain the eigenvalues and that we had previously, given by (27) and (28). The bottom right matrix is the same as , defined in (51), and so the eigenvalues of this matrix are and as defined in (52) and (53). The corresponding eigenvectors are given byWe now immediately get the linear stability of the steady state solution as all eigenvalues are strictly negative.

Lemma 5.7

The steady state (57) of Eqs. (20)–(23) is linearly stable when .

This result implies that all trajectories with initial conditions sufficiently close to the steady state will converge to it. On its own , it does not guarantee convergence for all initial conditions. However, for this model, a stronger stability result can again be proved.

Lemma 5.8

The steady state of Eqs. (20)–(23) is globally asymptotically stable when .

This result ensures that for any non-negative initial conditions and for all positive parameter values, each of the variables will converge to their unique steady state value, as with the previous model.

The proofs of these two results are again in Appendix A.

Rebound

In this section, we focus on the analysis of rebound in the latter stages of the evolution, i.e., the approach to the steady state. The linearised system will play an important role in this, but first we consider the limiting case and derive some estimates on the initial behaviour of the moderator.

For , the full TMDD equations with feedback via a moderator, i.e., (20)–(23), reduce toHence w is constant. With the initial condition (24), this implies and the system reduces to the TMDD equation without feedback as studied in our previous paper (Aston et al. 2014). However, if w has a different initial condition, say , then the system will converge to the steady state . Within the hyperplane , this equilibrium is an attractor. Thus if , then the receptor will converge to a value above the original baseline. Note that in this case, we have , corresponding to the invariance of the moderator w and the line of steady states .

For and the initial condition (24), we can show that initially the moderator w will increase to values above its baseline and hence the production of the receptor will be stimulated.

Lemma 5.9

For and , there exists a such that for all .

Proof

The initial condition is and since this is also the steady state value of w, we have . Differentiating (23) and evaluating at givesand the result follows immediately from this.

The initial conditions give that , hence initially y will decrease and take values smaller than its baseline value 1. The following lemma implies that if there is no rebound, then must be larger than 1 for all time. This observation will be used several times later to derive a contradiction on the assumption that there is no rebound.

Lemma 5.10

If , and for all , then for .

Assume that , and for all . We evaluate along the hyperplane givingBy hypothesis (h3), if and so in this case. Since for sufficiently small by Lemma 5.9, this result shows that cannot touch or cross the line when .

Next we will show by contradiction that is impossible too. We assume that there is some such that with for , i.e. is the first moment that equals one, and that also. Using the model equations (20)–(23), we obtainThus, using Taylor series, y and w near to are given byThe assumption for and for all together with imply that . Looking at the expressions for the derivatives at , this gives , and . Going to the next order in the Taylor series, since , then we must have that . If then this implies that for and sufficiently small, which contradicts our earlier assumption on w. Thus, the only possibility is that . Using the relation we then find thatand so this second derivative is zero if either or . If , then too. As we also have , this would imply that the solution is at the equilibrium at . Uniqueness of solutions gives that this is not possible. With the alternative option , evaluating , and on the lineshows that this line is invariant under the dynamics. At , , hence initially the solution is not on this line and therefore it can never be on the line in finite time, which rules out the possibility that and we have come to a contradiction. We therefore conclude that and so it follows that for all .

Corollary 5.11

If , and for some , then rebound occurs.

This result is simply the converse of Lemma 5.10.

The observation that initially the moderator will be larger than 1, combined with the fact that for , the initial condition leads to rebound, makes it likely that rebound will happen for all small values of . This is indeed the case as will be shown below.

Relative position of the eigenvalues

To analyse rebound in the latter stages of the evolution, we focus on solutions near to the baseline state. We are especially interested in the eigenvalue of the Jacobian matrix (58) with real part closest to zero, since generically trajectories approach the steady state tangent to the eigenvector corresponding to this eigenvalue. We note that the eigenvalues and are always real, depend on ,  and , do not depend on  and , and that . Furthermore, for and fixed, the eigenvalue is monotonic decreasing in  andOn the other hand, the eigenvalues and depend on and , do not depend on ,  and , and . Thus, the eigenvalue with real part closest to zero will be either or . As observed before, does not depend on , and for small we haveThus, for sufficiently small , the eigenvalue closest to zero will be . To describe how this changes when increases, we consider four intervals that are defined in terms of (see (56)) and (see Lemma 5.2)—a summary is sketched in Fig. 8: To visualise this information, we have sketched the intervals in the - plane in Fig. 8.

Fig. 8

The relative positioning of the eigenvalues in the – plane with and fixed. Recall that and , thus the three curves are relations of the form , with the appropriate function

(i.e., ):

Now is in the interval where takes real values, hence the eigenvalues and cross at the point where and solving this equation for gives Note that in this interval and so is well defined. Thus , when . Altogether, we conclude that in this interval, we have that if

(i.e., ):

Here is in the interval where takes complex values, hence for , is real and will be the eigenvalue closest to zero. At , and merge and become a complex conjugate pair. For just larger than , the real part of the complex pair will be closest to zero and this will be the case until the real part of the complex eigenvalues crosses . This occurs when and since , this crossing occurs at , where Since , it follows immediately that Note that , when . Thus in this interval, we have that for .

(i.e., ):

Now is in the interval where takes both real and complex values. If , then is real and we can use the calculation in the first point to conclude that if . If , then is complex and the calculation above gives that if . Note that in the open interval and that and collide with when . Furthermore, for . So altogether we have that in this interval

(i.e., ):

Since , we see immediately that in this interval for all .

Rebound in the latter stages of the evolution

Now that we have determined the positioning of the eigenvalues, we can state when rebound occurs in the latter stages of the evolution.

Theorem 5.12

Define m as in Theorem 4.10, i.e., and define , hence .These results are summarised in Fig. 9.

Fig. 9

Rebound overview, on the left for and on the right for . Note that . Recall that and

If , then rebound occurs generically for .

If and

, then rebound occurs generically for ;

, then rebound occurs generically for all .

If and

then rebound occurs generically for ;

then rebound occurs generically for .

The proof of this theorem will be given via two lemmas which determine when the eigenvalue closest to zero can be associated with rebound. First we consider the case that the eigenvalue (or its real part if complex) is closest to zero.

Lemma 5.13

If the eigenvalue (or its real part if complex) is closest to zero andthen rebound occurs (generically only for ).

or

, and ,

The proof of this lemma can be found in Appendix B.

Next we consider the case that the eigenvalue is closest to zero.

Lemma 5.14

If is the eigenvalue closest to zero, then rebound occurs generically if or if .

Again, the proof of this lemma can be found in Appendix B.

With the lemmas above, we are ready to prove Theorem 5.12.

Proof of Theorem 5.12

The results in this Theorem are essentially obtained by applying Lemmas 5.13 and 5.14 in each of the four regions described in Sect. 5.3.1. The results from this process are summarised in Table 1. The observations about the relative position of , and in the various regions follow from (61) and (62). In the regions , the case leads to the range . This range is not empty only when . The definition of shows that this condition holds if , which is why this condition is included in the second row. For the third row, we observe that , so the range is always empty if .

Table 1

Summary of results obtained by applying Lemmas 5.13, 5.14 to the eigenvalue categories described in Sect. 5.3.1. The ticks indicate that rebound occurs in the given intervals

Region	Re(λ_4)>λ_1	Conditions for rebound	Re(λ_4)<λ_1	Conditions for rebound
λ∗-≤λ_1<0	0<ϵ<ϵ_2(<ϵ1+)	✓	ϵ>ϵ_2	k_1>k_4 or ϵ_2<ϵ<-λ_1
λ_∞≤λ_1≤λ∗-	0<ϵ<ϵ_3(<ϵ1+)	✓	ϵ>ϵ_3	k_1>k_4 or ϵ_3<ϵ<-λ_1 if λ_1>-k_3
λ∗+≤λ_1<λ_∞	0<ϵ<ϵ_3(≤ϵ1+) ϵ>ϵ_2(≥ϵ1+)	✓ k_3>k_41+m	ϵ_3<ϵ<ϵ_2	k_1>k_4
λ_1<λ∗+	ϵ>0	0<ϵ<ϵ1+ or ϵ≥ϵ1+ if k_3>k_41+m	−	−

The two columns containing conditions for rebound in Table 1 can be combined in some cases, as shown in Table 2, where we separate out two regions in the parameter plane defined by and . For the third row, we have used that is monotonic decreasing in . Hence if , then (see (60)). Thus if and , then .

Table 2

Summary of conditions for rebound derived from Table 1

Region	k_1≤k_4	k_1>k_4
λ∗-≤λ_1<0	0<ϵ<ϵ_2, ϵ_2<ϵ<-λ_1	0<ϵ<ϵ_2, ϵ_2<ϵ
λ_∞≤λ_1≤λ∗-	0<ϵ<ϵ_3, ϵ_3<ϵ<-λ_1 if λ_1>-k_3	0<ϵ<ϵ_3, ϵ_3<ϵ
λ∗+≤λ_1<λ_∞	0<ϵ<ϵ_3	0<ϵ<ϵ_3, ϵ_3<ϵ<ϵ_2, ϵ>ϵ_2 if k_3>k_41+m
λ_1<λ∗+	0<ϵ<ϵ1+, ϵ≥ϵ1+ if k_3>k_41+m	0<ϵ<ϵ1+,ϵ≥ϵ1+ifk_3>k_41+m

We note in Table 2 that in many cases, there are neighbouring intervals with only the boundary point between them missing. Thus, filling in these specific points (which are proved below) gives the results shown in Table 3. The results in Theorem 5.12 are obtained by converting the eigenvalue ranges in Table 3 into conditions involving the parameters, using that and .

Table 3

Summary of final conditions for rebound

Region	k_1≤k_4	k_1>k_4
λ∗-≤λ_1<0	0<ϵ<-λ_1	ϵ>0
λ_∞≤λ_1≤λ∗-	0<ϵ<ϵ_3 or 0<ϵ<-λ_1 if λ_1>-k_3	ϵ>0
λ∗+≤λ_1<λ_∞	0<ϵ<ϵ_3	0<ϵ<ϵ_2 or ϵ>0 if k_3>k_41+m
λ_1<λ∗+	0<ϵ<ϵ1+ or ϵ>0 if k_3>k_41+m

We now consider a number of special cases that fill in the gaps that allow us to go from Tables 2 to 3.The rebound results that are summarised in Table 3 are shown in Fig. 9. There are two cases shown here, depending on whether or not the curve intersects the line . Since (see (60)), if then the curve is always below . Also, as (see (60)), then it can be shown that the line is above the line at this point, and so the line can never go below the curve as the curve is monotonically increasing with .

, .

When , then and so there are two repeated eigenvalues that are the least negative. With some calculations, it can be seen that the eigenvectors and collide (the denominator in the third and fourth entry vanishes, so the eigenvector has to be rescaled first to see this), hence we are in the case of an eigenvector and a generalised eigenvector. The eigenvector is and generically solutions will align with this eigenvector. In the proof of Lemma B.1, it is shown that when . Now for by (61) and this implies that the third and fourth components of have the same sign. The proof in Lemma B.1 now applies to prove that generically rebound occurs in this case.

, .

In this case there is triple eigenvector with one eigenvector together with two generalised eigenvectors. Again, generically solutions align with this eigenvector and, as above, the third and fourth entries have the same sign, thus and so rebound can be proved via a contradiction argument.

, with .

This is similar to the first case as the two eigenvalues and coincide when with an eigenvector and a generalised eigenvector in this case. Generically solutions will align asymptotically with the eigenvector. The proof of Lemma B.2 can now be applied in this case.

, with or .

In this case, the real part of the complex pair equals the eigenvalue . Approaching the steady state via the complex pair leads to rebound. The proof of Lemma 5.14 also applies with provided that and this occurs if and so rebound occurs if the steady state is approached tangent to . If the approach occurs via a linear combination of the three eigenvectors, then clearly rebound will then occur in this case also. Hence we can conclude that rebound will occur generically in this case.

Remarks

Some of the rebound results that we proved in Lemmas 5.13 and 5.14 hold only generically, and some cases always give rebound. However, for the sake of simplicity, we have not distinguished between these in Theorem 5.12 but simply refer to the existence of generic rebound.

In Theorem 5.12, we have proved the existence of rebound in a variety of cases. However, it should be noted that we have not proved or disproved the existence of rebound in other parameter regions, and so we do not have results for all possible parameter values.

From Theorem 5.12, we can conclude that there exists such that rebound occurs for all satisfying (see Fig. 9). Thus, for all parameter values, rebound occurs for sufficiently small . This might seem to contradict the result for the limit (the no-feedback case) where there is a large region in the parameter plane for which there is no rebound. The explanation is that for those parameter values, the magnitude of the rebound decays to zero as goes to zero.

Generalisations

Thus far, we have considered the basic TMDD model with feedback dynamics which is linear in the feedback moderator F. We now extend this model in two directions. First we extend the TMDD model to allow for more compartments and secondly we consider more general types of feedback dynamics. Below we show that some of the rebound results can be extended to these more general models.

Generalised TMDD model with feedback via a moderator

In this section we consider a generalisation of the basic TMDD model to more compartments. It is assumed that the receptor is still only in one compartment. In the basic TMDD model, the variable y represents the normalised receptor in this compartment. Now we will assign y to represent either the normalised free receptor (as before) or it could be the normalised total amount of receptor. The normalised ligand-product components (x and z) of the basic TMDD model are generalised to a vector that represents the ligand and product in the compartments. The generalised TMDD model with feedback becomes where , and . We note that if depends only on but not on y, or if depends only on y but not on , then the coupling between and y is in one direction only. In this case, the model could not be classed as a generalised TMDD model, although our results below still hold in these cases. We assume that the assumptions (h1)–(h4) on the feedback function h(y) hold and we make the following assumptions on the functions and :With these assumptions, the Jacobian matrix evaluated at the steady state is given byNote that the condition leads to a block lower triangular structure in this matrix. This implies that the eigenvalues of this matrix are the eigenvalues of the Jacobian matrix and the eigenvalues of the bottom right matrix in . The global asymptotic stability of the steady state implies that all eigenvalues have negative or zero real part. Note that the eigenvalues of do not depend on  and that our final assumption implies that the eigenvalue closest to zero will have strictly negative real part. The bottom right matrix is the same as in the basic TMDD model with feedback. We denote its eigenvalues by and , where the eigenvalue corresponds to earlier and to . Thus is the closest eigenvalue to zero for sufficiently small .

There exists an such that the Eqs. (63)–(65) have the steady state solution , . This implies that and . We also assume that this steady state is globally asymptotically stable.

In the basic TMMD model, a change in the receptor level while the ligand and product are at baseline does not move the ligand or product away from baseline and we make a similar assumption for this model, i.e., for all which implies that . Similarly, the function is not affected by the receptor dynamics if the ligand and product are at baseline, i.e., for all which implies that . In fact, in both cases, all that we really require is the weaker derivative condition.

The initial conditions are given by for some with such that . Hence initially the amount of receptor will decrease.

The Jacobian matrix has no zero or purely imaginary eigenvalues.

We definewhere is the value of at which the eigenvalue coincides with the eigenvalue closest to zero of the matrix . Recall that is the value of for which and become a complex pair of eigenvalues. We then have the following result.

Theorem 6.1

If the above assumptions on and hold, then rebound occurs generically for .

The condition ensures that , which is required in the proof of Lemma 5.9. The w equation is unchanged in this setting and so Lemma 5.10 also holds. The proof of Lemma B.1 also holds in this setting, although the exceptional case that the trajectory approaches the steady state tangent to a different eigenvector cannot be excluded outright in this more general case, and so we only have a generic result.

If the bottom right matrix in has complex eigenvalues with real part closest to zero of all the eigenvalues, then again there will be oscillation about the steady state in y as the steady state is approached and so rebound occurs infinitely often but with exponentially decreasing amplitude.

Generalisation of the dynamics of the feedback

In our considerations so far, the feedback dynamics has been linear in the feedback moderator F (or w in the non-dimensional equations), see (14) or (23). In this section we will consider a more general form of feedback, and focus on the case of a small value of the parameter (). In particular, we again assume that the receptor equation is described by (13) and describe the dynamics of F byNon-dimensionalising as before, we obtain a new equation for the feedback variable w () given bywhereand is defined by (25). We make a number of assumptions regarding the function g:The assumption (h1’) ensures that the equations retain the steady state with . Assumption (h2’) ensures that when the system is perturbed from the steady state by reducing y, which is what happens when the ligand is injected, then w will increase, resulting in greater production of the receptor. Finally, assumption (h3’) ensures that the steady state of (67), assuming that is held fixed, is stable.

We now show that some of the results obtained previously also hold in this more general setting. We again start by considering the Jacobian matrix of Eqs. (20)–(22), (67) evaluated at the steady state (57), which is given by

Theorem 6.2

If assumptions , and hold, and if for all , then there exists such that rebound occurs in the model given by Eqs. (20)–(22), (67) for all satisfying .

We note that assumptions (h1) and (h2) imply that for sufficiently close to 1. The extra assumption requires that this holds for all .

The proofs of Lemmas 5.9, 5.10 and B.1 can easily be adapted to this modified model. In particular, for Lemma 5.9, as before and in this case, we haveand so Lemma 5.9 holds.

For Lemma 5.10, we note thatusing our stated assumption. Combined with Lemma 5.9, this implies that cannot touch or cross the line when . The second part of the proof of Lemma 5.10 shows by contradiction that is impossible and the same proof holds in this case with replaced by .

The Jacobian matrix associated with Eqs. (49) and (67) (analogous to (51)) is given bywhich has eigenvalues , when as previously. These two eigenvalues collide when the discriminant of the characteristic polynomial is zero, and this first occurs whenWith this definition of , we claim that Lemma B.1 holds also. To see this, we note that the eigenvector given by (59) is the same in this case, but with as the eigenvalue of (69) that is closest to zero. We now haveIt is easily verified that and this implies that for . Thus, if , then the third and fourth components of the eigenvector have the same sign and the proof of Lemma B.1 again implies that rebound occurs.

The second part of the proof of Lemma B.1 considered the non-generic cases of the trajectory approaching the steady state tangent to one of the other eigenvectors. For the eigenvector in (59), we note that it is easily verified that and so the same arguments apply in this case. The argument for the eigenvector is unchanged, while for , the requirement for the third and fourth components of the eigenvector to be the same is and it is easily verified that this holds.

If corresponds to the point at which the eigenvalues and coincide, then the conditions of Lemma B.1 hold with and so we conclude that rebound therefore occurs for all satisfying .

This result shows that for this generalised model, rebound always occurs for sufficiently small even if there is no rebound without feedback (), as is the case for the original TMDD model with feedback moderator.

Clearly, if the bottom right matrix in (68) has complex eigenvalues which have real part closest to zero of all the eigenvalues, then there will again be oscillation about the steady state in y and so rebound occurs infinitely often but with exponentially decreasing magnitude.

Applications

We now consider two applications of these results to two different models. The first is the standard TMDD model with feedback. The second is a more complicated model of the effect of efalizumab on patients with psoriasis. This model is an application of the generalizations in Sect. 6.

Example 1

We consider the TMDD equations including feedback with moderator given by (11)–(14). We use parameter values for the IgE mAb omalizumab (Sun 2001; Agoram et al. 2007), namely , , , , , , , and initial drug dose  nM. These parameter values have and so no rebound will occur when there is no feedback or in the case of direct feedback for any feedback function H(R). This follows directly from Theorem 2.1, Corollary 4.4, and Theorem 4.10 since implies that the parameter values are in the left green region in Figs. 3, 4, and 6.

However if there is feedback via a moderator then Theorem 5.12 implies that rebound occurs for all sufficiently small , even though there is no rebound at . The parameters for this example giveThus Theorem 5.12 gives that the parameters are in the upper left region of Fig. 9 and that rebound will happen for (using ). However, this is a lower bound on the region of rebound. To illustrate this, we have used the mainly linear feedback function H(R) in (34) with . The magnitude of the rebound is sufficiently small in this case that a nonlinear function is not required. The maximum rebound is plotted in Fig. 10 for a range of values of which does indeed show that there is rebound for all and sufficiently small. For this example, the rebound ends at approximately . In this case, a local maximum for R persists for larger values of , but with . The value of time at which the maximum rebound occurs is also shown in Fig. 10. We note that as then and so the time of maximum rebound moves to infinity as the magnitude decreases. The same mechanism does not apply when the rebound disappears around since a local maximum persists in this case but at values lower than baseline, as already discussed.

Fig. 10

Left a plot of the relative magnitude of rebound for feedback with a moderator and the mainly linear feedback function (34) with for a range of values of . Right a plot of the time at which the maximum rebound occurs

We have numerically sampled larger values of and there is no indication of further rebound occurring. Of course we know that when , the feedback with moderator becomes the direct feedback case, for which there is no rebound in this example.

Example 2

In a study on the effect of efalizumab on patients with psoriasis (Ng et al. 2005), rebound was reported in some patients. The model equations proposed in Ng et al. (2005) arewhere , and are the amounts of efalizumab in the depot, central and peripheral compartments respectively, is the total %CD11a on the surface of each T cell and is the production rate of %CD11a to the T cell surface. In this section we will show that this model can be written as a generalised model with feedback as studied in Sect. 6.1. Strictly speaking, this model is not a TMDD model as the efalizumab equations decouple from the %CD11a equations. As indicated in Ng et al. (2005), one could use instead of in the PK equations and have a more TMDD-like model. In Ng et al. (2005), the simpler model was chosen as there was no difference in the fit of the data. If one would use the more complicated model, the parameters would need to be fitted again. So we focus in this section on the same model as in Ng et al. (2005). However, our analysis also applies to the more complicated model and the observed rebound would again be due to the slow feedback response.

In this model, is the total receptor, free and bound, and is the feedback moderator. The plots in Fig. 3 of Ng et al. (2005) show rebound in the total receptor, which of course is not the same as rebound in the free receptor. However, we first show that this model fits with the generalised model that we considered in Sect. 6.1 and we will then discuss the issue of rebound in the free receptor. To put this model into our framework, we take , and . The differential equations for Y and F are then given by whereThe non-dimensional form of these equations is precisely in the form of Eqs. (63)–(65). These equations have the unique steady state in the region of interest given bywhere is the positive root of the quadratic equationSubstituting into (72) gives a quadratic equation for which has two negative roots, and this implies that . It is then easily verified that all the conditions stated in Sect. 6.1 on the functions and hold. The particular function H that is used in this model is given byWe non-dimensionalise, considering only Y and F, by defining the new variablesThe non-dimensional form of equations (70)–(71) is thenwhere dot denotes differentiation with respect to , is the (unspecified) non-dimensional form of , andNote that we have used (72) in deriving h(y). All of our assumptions (h1)–(h4) on the function h(y) are now satisfied andsince . We tried to reproduce the numerical results shown in Fig. 3B (dashed line) of Ng et al. (2005) but found that this was only possible by changing from the stated value of 0.00154 day to 0.0154 day (the smaller value of did show rebound, but with a smaller magnitude of about 110% of baseline). Moreover, the authors in Ng et al. (2005) report an affinity of 0.033 g/ml, which is M = M, or 200pM. If we use the reported of 0.00154 day and assume a typical (see Aston et al. 2011) of 0.592 nM/day, the ratio pM, i.e. 100-fold more potent than reported. A ten-times higher still does not give the right affinity but at least it is closer. The values of the other parameters used, from Tables II and III of Ng et al. (2005), are  day,  day,  day,  day,  day,  mL/kg, g/mL,  day, g/mL, ,  %CD11a/day. The plot of the total receptor in this case is the solid blue line shown on the left in Fig. 11.

Fig. 11

The total and free %CD11a relative to baseline after a single 3 mg/kg intravenous dose of efalizumab. On the left is the plot with feedback turned on and rebound can be observed. On the right is the plot without feedback and no rebound occurs

To find out how feedback influences the model, we ran a simulation of the model when the feedback was turned off, i.e., we took to be constant at its baseline value. As can be seen on the right in Fig. 11, without feedback there is no rebound. So clearly the rebound in this example is caused by the feedback.

Next we will interpret these results using our analysis. Using the parameter values in Ng et al. (2005), but the modified value of , we find thatand since this is so small, we expect that the eigenvalue closest to zero of the Jacobian matrix evaluated at the steady state will come from the lower diagonal block. This is indeed the case since the eigenvalue closest to zero of this block is while the eigenvalue closest to zero from the remaining block is . Thus, in the notation of Theorem 6.1, we have and so this Theorem says that we generically have rebound in the total receptor , as confirmed in Fig. 11 and by the numerical results in Ng et al. (2005).

The two eigenvalues of the lower diagonal block collide at when , but since there are eigenvalues of that are closer to zero than this, we would not expect to see oscillatory convergence to the steady state in this model by simply increasing (), as described in Lemma B.3.

Finally, we recall that the variable in the above model is the total of the free and bound receptor. The free receptor on T cells is given by (Ng et al. 2005) and we note that if , then the free receptor is the same as the total receptor. Now is the amount of efalizumab in the central compartment and while this is present, the total receptor is kept at a low level. However, once has been almost depleted, then starts to rise again and rebound occurs. Thus, is very low in the phase of the dynamics in which rebound in occurs, and this implies that the free receptor is very similar to and so will also rebound. The free receptor for the above example is also shown in Fig. 11 as the red dashed line and it can clearly be seen that this also has rebound which is very similar to that for the total receptor as anticipated.

Conclusions and discussion

In this paper we have extended our previous work on rebound in the basic TMDD model by including feedback. The feedback is negative, thus the rate of receptor production increases when the concentration drops below baseline (and vice versa). We modelled the feedback with an additional dynamic equation for the feedback moderator. If the feedback responds very fast, a quasi equilibrium approximation can be used and the feedback can be incorporated into the synthesis term itself.

We have shown that rebound is more likely to occur if feedback is present and that the likelihood of feedback occurring depends on the response speed of the feedback. If the feedback responds sufficiently slowly there will always be rebound. In the psoriasis example presented in Sect. 7.2, it can be seen that the rebound can be significant, an increase of over 40% of baseline is observed. The TMDD example in Sect. 7.1 shows that this is not always the case. In general, it is a challenge to describe the main parameters that influence the magnitude of the rebound. One important parameter is obviously the response speed of the feedback. If the response of the feedback goes to zero, then we converge to the case without feedback. Thus for many values of the elimination parameters the magnitude of the rebound will diminish if the response of the feedback slows down to zero and in the limit the rebound will have gone. The exceptions are of course those elimination parameter values for which rebound is observed if no feedback is present (the elimination rate of the antibody-protein complex is less than both the elimination rate of the antibody and the elimination rate of the protein ( and )).

At first sight, it might be surprising that there is always rebound for a very slow feedback response rate (), while there are large regions without rebound when the feedback moderator does not change and there is no feedback (). However, it can be seen that, in the absence of feedback dynamics (), a larger initial condition for the moderator w will lead to a steady state with which is above the baseline. Furthermore, Lemma 5.9 shows that when the moderator starts at baseline, ligand is added and feedback is turned on (), then the moderator w initially increases to above its baseline value, albeit on a slow time scale since is small. However, this implies that w will also only decrease back to its steady state value of 1 on a slow time scale. While , y converges towards w on a faster time scale and will therefore itself exceed 1, resulting in rebound, followed by a slow relaxation of y and w back to their final steady state.

For moderately slow feedback response rates, the rebound is expected to be maximal, while for very fast feedback responses (), the rebound will be approximated by the direct feedback limit. The results for the direct feedback model are qualitatively very similar to the results without feedback: the existence or non-existence of rebound again depends only on the three elimination rates , and . An extra feature now is that the characteristics of the feedback function H(R) play a role as well if . In particular, in both the direct feedback and no feedback models, there is no rebound if the elimination rate of the ligand is less than the elimination rate of the antibody-protein complex (). Furthermore, in both models there will be rebound if the elimination rate of the antibody-protein complex is less than both the elimination rate of the antibody and the elimination rate of the protein ( and ). However, if , then the presence of feedback will increase the region in which rebound occurs. The details depend on the characteristics of the feedback function H(R) and can be found in Theorem 4.10. Again, in both models, the association and dissociation rate parameters and do not influence the occurrence of rebound, though they will play a role in the magnitude of the rebound if it occurs. There are some preliminary observations that a large will dampen rebound, but this needs further work.

In the general feedback model, we have focused on proving the existence of rebound. In general, it is challenging to prove the non-existence of rebound as it has to be shown that the receptor is bounded by its baseline for all time. We expect that no rebound will occur for or small for a sufficiently fast response rate ( large), but we have not yet obtained analytical results about this. The fast response is a singular limit and it can be shown that the dynamics has to be close to the limiting case. Thus if any rebound did happen, it would have to be very small.

Furthermore, we have also shown in this paper that most of the results obtained on rebound in the model with slow feedback can be extended to more general systems of TMDD equations such as equations with more compartments. We illustrated this with the example of psoriasis, where rebound is predicted by the model and underpins the observations in patients (see Sect. 7).

It must be acknowledged that during target discovery and validation, knowledge of the quantitative relationship between the target neutralisation and the nature of the negative feedback is often not known. However, using prior knowledge on the target of interest and similar targets, the model can be used to evaluate the possibility of feedback and its potential magnitude, and thus confirm the suitability of the target.