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.
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 On the horizontal axis is the elimination rate of the ligand (antibody/drug), denoted by k e(L) , and on the vertical axis is the elimination rate of the receptor (antigen/target), denoted by k out . The elimination rate of the antibody-antigen complex (drugtarget product) is denoted by k e(P) . 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, m = h = M = 0 and the region in which there is rebound is the red region above the line k out = k e(P) . 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)

Fig. 2
The TMDD reaction mechanism 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 k on and k off , the receptor turnover and elimination rates k in and k out , and the elimination rates of the ligand and complex k e(L) and k e(P) . The system is assumed to be initially at steady state, into which a single bolus infusion L 0 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 L = P = 0, R = k in /k out . Adding the bolus injection L 0 gives the initial conditions After 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 R 0 . In Aston et al. (2014), we determined precise conditions for the existence and non-existence of rebound. Our main result was the following.

Theorem 2.1 Rebound occurs in Eqs. (1)-(3) if and only if
k e(P) < k e(L) and k e(P) < k out .
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 (k e(L) , k out ) plane, as shown in Fig. 3. The above condition for rebound in terms of non-dimensional parameters is given by k 4 < k 1 and k 4 < k 3 , Fig. 3 The green region shows the part of the (k e(L) , k out ) parameter plane where there is no rebound, and the red region shows where there is rebound. The boundaries are given by the lines k e(L) = k e(P) and k out = k e(P) . This figure is similar to Fig. 2 in Aston et al. (2014)  where 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 (k in ) is not constant but varies as R changes. Thus, Eq. (2) is now modified to have the form The 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 equation The constant α ≥ 0 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 α = 0, 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 relation is 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: • the production rate of the receptor, denoted by k in F, should be positive; • the receptor concentration should have the baseline R = R 0 = k in /k out ; • 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.
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 F > 0, 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 that d F dt F=0 > 0 and so the line F = 0 cannot be crossed from above to below. This also ensures that F > 0 for all t > 0 and hence the production rate is always positive. The basic TMMD model without feedback, Eqs.
(1)-(3), has a steady state solution given by L = P = 0, R = R 0 = k in /k out . 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 F = 1. 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 H (R 0 ) ≤ 0. We note that if H is a positive, strictly decreasing, smooth function with H (R 0 ) = 1, 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 F(0) = 1 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 variables We also define a new function The above assumptions on the function H (R) translate into the following assumptions on h(y): (h1) h(y) > 0 for all y ≥ 0; (h2) h(1) = 1; (h3) h(y) > 1 when 0 ≤ y < 1 and 0 < h(y) < 1 when y > 1; and note that the assumptions (h2)-(h4) imply that h 0 ≥ 0. In terms of the non-dimensional quantities, the direct feedback equations (8)-(10) becomeẋ with initial conditions where dot denotes differentiation with respect to τ and the dimensionless parameters are defined as in Aston et al. (2014) by Clearly, this choice of non-dimensionalisation requires that k on = 0 and k in = 0 (so that R 0 = 0). 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 conditions Furthermore, = α 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 x, y, z ≥ 0 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. 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 R 3 defined by x, y, z ≥ 0, the Eqs. (16)-(18) have a unique steady state given by This 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 λ 1 and λ 2 are the same as previously and are given by while the third eigenvalue, which previously was λ 3 = −k 3 , is now given by where we recall that h 0 = −h (1). 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 x = z = 0 satisfies (16) and (17). Substituting x = z = 0 into (18), we obtain the single differential equatioṅ Assumption (h2) ensures that the previous steady state value y = 1 is also a steady state of this equation while assumption (h3) ensures that this is the only non-negative steady state solution and thatẏ < 0 if y > 1 andẏ > 0 if 0 ≤ y < 1. The theory of scalar, autonomous differential equations (Jordan and Smith 2007) then tells us that if the initial condition y 0 is less than the steady state (y 0 < 1) then y(τ ) will increase monotonically to the steady state, whereas if y 0 is greater than the steady state (y 0 > 1), then y(τ ) will decrease monotonically to the steady state. As y(τ ) 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 that As h 0 ≥ 0, the inclusion of the function h(y) ensures that the asymptotic rate of convergence to the steady state (proportional to e −k 3 (h 0 +1)t ) will be faster than or equal to the rate of convergence which occurs when h(y) = 1 (proportional to e −k 3 t ). 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 y 0 < 1, then y(τ ) < 1 for all τ ≥ 0.

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 y ≤ 1. 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 h 0 > 0, 0 < β < 1, and This function satisfies all the conditions (h1)-(h4). For the linear part of the function h, we note that h(y) − y = (1 + h 0 )(1 − y) and so Eq. (18) becomeṡ We then have the following result. Proof We note that for 0 ≤ y ≤ 1, Eqs. (16), (17) and (33) are the same as the basic TMDD model that we studied in Aston et al. (2014) except that k 3 has been replaced by k 3 (1 + h 0 ). 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 0 ≤ y ≤ 1 and so the matching point y = 1 + β/ h 0 in (33) is not reached and the feedback function for y > 1 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 (x, y, z) = (0, 1, 0) and the dynamics for 0 ≤ y ≤ 1 + β h 0 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 k 3 replaced by Converting back to the dimensional parameters gives the following result.

Corollary 4.4 Consider the mainly linear feedback function
, H 1 (R) ∈ (0, 1) and H 1 and its first and second derivatives match the linear function at R = R 0 + β/H 0 .
Rebound occurs in the model given by Eqs. 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.

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 form Differentiating (35), evaluating at y = 1 and using (15), we see that Condition (h3) implies that h(y) > 0, when 0 ≤ y < 1, and 0 <h(y) < 1 y − 1 , when y > 1.
With this form of h, Eq. (18) becomeṡ Remark We note that expanding h(y) in a Taylor series with remainder gives that h(y) = −h (ξ(y)) where ξ(y) lies between y and 1. However, the fact that ξ(y) 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. Proof Following through the proof of Theorem 3.4 of Aston et al. (2014) for these modified equations, the only difference is thaṫ Nowh(0) > 0 by (37) and soẏ| y=0 > 0 when z ≥ 0, 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
Proof 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 k 3 replaced by k 3 (1 + h 0 ). 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 k 3 . Since λ 1 does not involve k 3 , there is no change to this eigenvalue.
The proofs in regions I and IV of the (k 1 , k 3 ) parameter plane in Aston et al. (2014) both made use of the fact that for the model without feedback the plane v = y + z = 1 is invariant when k 3 = k 4 . 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 v = y + z is now given bẏ and sov Clearly no condition on the parameters will give this derivative to be zero (unless h is constant, the case considered in the previous section), and so we see, as already stated, that the plane v = 1 is no longer invariant. However, the function h(y) has a continuous derivative by (h4) and this implies thath(y) is continuous. Thus the Extreme Value Theorem for continuous functions gives thath(y) has a minimum and a maximum for y ∈ [0, 1]. So we define Using (36), this immediately implies that Also, sinceh(y) > 0 for all 0 ≤ y < 1 by (37), we get Since m and M are the minimum and maximum ofh(y), then clearly m ≤h(y) ≤ M for all 0 ≤ y ≤ 1, and substituting forh(y) from (35) and rearranging gives Graphically, this means that on the interval y ∈ [0, 1], the function h(y) is bounded by two straight lines, both of which pass through the point (1, 1), and which have slopes −m and −M. This situation is illustrated in Fig. 5. Note that the Mean Value Theorem gives that −M ≥ min{h (y) | 0 ≤ y ≤ 1} and −m ≤ max{h (y) | 0 ≤ y ≤ 1}. 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 (k 3 , k 4 )-plane for which the dimensionless variable v associated with the total amount of receptor stays on one side of the plane v = 1.
Proof We consider the vector field at the plane v = 1. Equation (40) giveṡ Since we have restricted attention to the plane v = 1, this implies that y + z = 1 and so y ≤ 1 as z must be non-negative. As 1 − y is non-negative, we can use the bounds onh(y) given in(41) to get First we consider the case , hence the vector field on the plane v = 1 is either tangent to v = 1 or points towards v < 1. This implies that the region of phase space defined by x, y, z ≥ 0, v = y + z ≤ 1 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 x = 0, y = 0 and z = 0 in Lemma 4.1 and for the plane v = 1 above. Since v(0) = 1 and this region is invariant, we conclude that v(τ ) ≤ 1 for all τ ≥ 0. Next we consider the case k 4 ≤ k 3 (1 + m). Now the other inequality in (45) giveṡ With this observation, we can now derive a result on rebound for regions similar to regions I and IV in Aston et al. (2014). (41), i.e., the function h(y) satisfies (44).

(i) Rebound does not occur if
Using (42) and (43), we note that This 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 (i) The case when (46) holds is similar to region I in Aston et al. (2014) and has the property that v(τ ) ≤ 1 for all τ ≥ 0 by Lemma 4.7. Thus for all τ ≥ 0, we have y(τ ) ≤ v(τ ) ≤ 1 and so rebound cannot occur in this case. (ii) When (47) holds (a condition similar to the one for region IV in Aston et al. (2014)), Lemma 4.7 ensures that v(τ ) ≥ 1 for all τ ≥ 0. The right hand inequality in (47) implies that λ 1 < λ 3 = −k 3 (1 + h 0 ) and so the eigenvalue closest to zero is λ 3 , 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 λ 3 points along this axis. As the trajectory always satisfies v ≥ 1 for all τ ≥ 0 then we again conclude that generically, almost all orbits will approach y = 1 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 v 2 is outside of the phase space and so an orbit cannot approach tangent to this onedimensional manifold.
The one-dimensional linearised manifold in the direction of the eigenvector v 1 given in (30) is We note that λ 1 + k 3 (1 + h 0 ) < 0 from the right hand inequality in (47). Thus, we must take a < 0 to ensure that this linearised manifold gives positive values of x and z (since k 2 + k 4 + λ 1 > 0 by Lemma 2.7 of Aston et al. 2014). The linearised manifold for v = y + z 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 v < 1. Since the trajectory in this case must satisfy v ≥ 1, 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 λ 1 = λ 3 with k 3 > k 4 /(1 + m). We now show that there is also rebound along this line.

Theorem 4.9 Rebound occurs if k
Proof The proof uses exactly the same arguments as the proof of Theorem 3.7 in Aston et al. (2014). Along the curve k 3 = −λ 1 /(1 + h 0 ), the eigenvalues λ 1 and λ 3 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 k 4 /(1 + m) < k 3 , Lemma 4.7 gives that v(τ ) ≥ 1 and the intersection of the y-axis and the region v(τ ) ≥ 1 is the part of the y-axis with y ≥ 1. Hence all orbits must approach the steady state (x, y, z) = (0, 1, 0) from above and rebound will occur.

For the model equations (16)-(18) and for any function h(y) satisfying the assumptions
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 define These results are illustrated in Fig. 6. We now focus on a few aspects of the theorem and compare it with the case where there is no feedback (F = 1, and so h 0 = 0 = H 0 ).
• In the equations without feedback, rebound occurred in the region k e(P) < k e (L) and k e(P) < k out (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 k e(L) ≤ k e(P) (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 k out < k e(P) < k e(L) (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.
we have no general results, and the existence or otherwise of rebound will depend on the particular function h.
In this case, the region of uncertainly disappears and we get that rebound will occur if and only if k e(P) < k e(L) and k e(P) < k out (1 + H 0 R 0 ). • If H is nonlinear but with m = R 0 H 0 (in which case the lower bounding curve for H (R) is the tangent to the curve at R = R 0 ), then the red region where there is rebound for k out > k e(P) /(1 + R 0 H 0 ) in Fig. 4 is preserved and the region of uncertainty is below this line. Conversely, if M = R 0 H 0 (so that the upper bounding curve for H (R) is the tangent at R = R 0 ), then the green region where there is no rebound for k out ≤ k e(P) /(1 Fig. 4 is preserved and the region of uncertainty is above this line. • We note that when k e(L) < k e(P) (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 k e(L) > k e(P) , then the form of the feedback function determines whether or not rebound occurs.
In summary, the region where rebound occurs increases slightly to include some more values with k out < k e(P) < k e(L) when direct feedback is introduced into the basic TMDD model.

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 x = z = 0 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 y = w = 1, i.e. the baseline values. Below we will prove the following rebound result.
Theorem 5.1 For any ∈ (0, ∞), there are always some initial conditions with y(0) < 1 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 . The eigenvalues of J 0 can be found as solutions of the quadratic characteristic polynomial derived from this matrix. The discriminant of this quadratic equation is and the eigenvalues are while the eigenvectors are of the form Note that in this section, the definition of the eigenvalue λ 3 is different from that used in the previous section. The discriminant D is itself a quadratic function of which Movement of the eigenvalues λ 3 and λ 4 in the complex plane for increasing (indicated by the direction of the arrows) has two positive roots given by Clearly the discriminant is negative between these roots and positive elsewhere, and so we have the following result.
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 .
When = 0, we see that λ 3 = −k 3 and λ 4 = 0. As increases from zero, the two eigenvalues move towards each other and collide on the negative real axis when = − 1 . After this collision, the eigenvalues become complex with Re(λ 3,4 ) = −(k 3 + )/2, which decreases monotonically with . As increases further, the two complex eigenvalues collide on the real axis when = + 1 and become real again. As increases from + 1 to ∞, λ 4 increases from the point of collision to −k 3 (1 + h 0 ) and λ 3 decreases to −∞ from this point. At the collision values = ± 1 , the eigenvalues λ 3 = λ 4 are given by Note that this implies that λ + * < λ ∞ < λ − * , as λ ∞ = −k 3 (1 + h 0 ). 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.
Thus the trajectory for all non-negative initial conditions will eventually approach the steady state (y, w) = (1, 1) 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.
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.

Invariance, steady states and stability
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 x, y, z, w ≥ 0 for Eqs. 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 R 4 defined by x, y, z, w ≥ 0, 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 ⎛ (1)). The eigenvalues of this matrix can be found from the top left 2 × 2 matrix and the bottom right 2 × 2 matrix. From the top left matrix, we obtain the eigenvalues λ 1 and λ 2 that we had previously, given by (27) and (28). The bottom right matrix is the same as J 0 , defined in (51), and so the eigenvalues of this matrix are λ 3 and λ 4 as defined in (52) and (53). The corresponding eigenvectors are given by (59) We now immediately get the linear stability of the steady state solution as all eigenvalues λ i are strictly negative.
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
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 = 0 and derive some estimates on the initial behaviour of the moderator. For = 0, the full TMDD equations with feedback via a moderator, i.e., (20)-(23), reduce toẋ Hence w is constant. With the initial condition (24), this implies w = 1 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 w = w 0 = 1, then the system will converge to the steady state (x, z, y, w) = (0, 0, w 0 , w 0 ). Within the hyperplane w = w 0 , this equilibrium is an attractor. Thus if w 0 > 1, then the receptor will converge to a value above the original baseline. Note that in this case, we have λ 4 = 0, corresponding to the invariance of the moderator w and the line of steady states (0, 0, w 0 , w 0 ). For > 0 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.
Proof The initial condition is w(0) = 1 and since this is also the steady state value of w, we haveẇ(0) = 0. Differentiating (23) and evaluating at τ = 0 gives and the result follows immediately from this.
The initial conditions give thatẏ(0) = − 1 μ , 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 w(τ ) 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.
Going to the next order in the Taylor series, since y(τ ) ≤ 1, then we must have thatÿ(τ ) ≤ 0. Ifÿ(τ ) < 0 then this implies that w(τ +τ ) < 1 forτ < 0 and |τ | sufficiently small, which contradicts our earlier assumption on w. Thus, the only possibility is thatÿ(τ ) = 0. Using the relation x(τ ) = μk 2 z(τ ) we then find thaẗ and so this second derivative is zero if either z(τ ) = 0 or k 1 = k 4 . If z(τ ) = 0, then x(τ ) = μk 2 z(τ ) = 0 too. As we also have y(τ ) = 1 = w(τ ), this would imply that the solution is at the equilibrium at t =τ . Uniqueness of solutions gives that this is not possible. With the alternative option k 1 = k 4 , evaluatingẋ − μk 2ż ,ẏ andẇ on the line shows that this line is invariant under the dynamics. At t = 0, x(0) − μk 2 z(0) = 1, 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 k 1 = k 4 and we have come to a contradiction. We therefore conclude that w(τ ) = 1 and so it follows that w(τ ) > 1 for all τ > 0.
Proof 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 = 0, the initial condition w(0) > 1 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 f o r a ll ε > 0 Fig. 8 The relative positioning of the eigenvalues in the k 1 -k 3 plane with k 2 and k 4 fixed. Recall that , thus the three curves are relations of the form k 3 = λ 1 (k 1 , k 2 , k 4 )/ f (h 0 ), with f (h 0 ) the appropriate function the eigenvalues λ 1 and λ 2 are always real, depend on k 1 , k 2 and k 4 , do not depend on and k 3 , and that λ 2 < λ 1 . Furthermore, for k 2 and k 4 fixed, the eigenvalue λ 1 (k 1 , k 2 , k 4 ) is monotonic decreasing in k 1 and On the other hand, the eigenvalues λ 3 and λ 4 depend on k 3 and , do not depend on k 1 , k 2 and k 4 , and Re(λ 3 ) ≤ Re(λ 4 ). Thus, the eigenvalue with real part closest to zero will be either λ 1 or λ 4 . As observed before, λ 1 does not depend on , and for small we have Thus, for sufficiently small , the eigenvalue closest to zero will be λ 4 . To describe how this changes when increases, we consider four λ 1 intervals that are defined in terms of λ ± * (see (56)) and λ ∞ (see Lemma 5.2)-a summary is sketched in Fig. 8: ): Now λ 1 is in the interval where λ 4 takes real values, hence the eigenvalues λ 4 and λ 1 cross at the point where λ 4 = λ 1 and solving this equation for gives Note that λ 1 + k 3 (1 + h 0 ) = λ 1 − λ ∞ = 0 in this interval and so 2 (λ 1 , k 3 ) is well defined. Thus 2 (λ 1 , k 3 ) = − 1 (k 3 ), when λ 1 = λ − * . Altogether, we conclude that in this λ 1 interval, we have that Here λ 1 is in the interval where λ 4 takes complex values, hence for < − 1 (k 3 ), λ 4 is real and will be the eigenvalue closest to zero. At = − 1 (k 3 ), λ 4 and λ 3 merge and become a complex conjugate pair. For just larger than − 1 , 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 λ 1 . This occurs when Re(λ 3,4 ) = λ 1 and since Re(λ 3,4 ) = −(k 3 + )/2, this crossing occurs at = 3 , where Note , then λ 4 is complex and the calculation above gives that Re(λ 4 ) > λ 1 if 0 < < 3 (λ 1 , k 3 ). Note that in the open λ 1 interval and that 3 (λ 1 , k 3 ) and 2 (λ 1 , k 3 ) collide with + 1 (k 3 ) when λ 1 = λ + * (k 3 ). Furthermore, 2 (λ 1 , k 3 ) → ∞ for λ 1 → λ ∞ . So altogether we have that in this interval , we see immediately that in this interval Re(λ 4 ) > λ 1 for all > 0.
To visualise this information, we have sketched the λ 1 intervals in the k 1 -k 3 plane in Fig. 8.

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. and on the right for k 2 > k 4 .
These results are summarised in Fig. 9.
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 λ 4 (or its real part if complex) is closest to zero.

Lemma 5.13
If the eigenvalue λ 4 (or its real part if complex) is closest to zero and • 0 < < + 1 or • ≥ + 1 , and k 3 > k 4 1+m , then rebound occurs (generically only for − 1 < < + 1 ). The proof of this lemma can be found in Appendix B.
Next we consider the case that the eigenvalue λ 1 is closest to zero.
Lemma 5.14 If λ 1 is the eigenvalue closest to zero, then rebound occurs generically if k 1 > k 4 or if < −λ 1 .
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  Table 2 Summary of conditions for rebound derived from Table 1 Region results from this process are summarised in Table 1. The observations about the relative position of 2 , 3 and + in the various regions follow from (61) and (62). In the regions λ + * < λ 1 < λ − * , the case Re(λ 4 ) < λ 1 leads to the range 3 < < −λ 1 . This range is not empty only when 3 < −λ 1 . The definition of 3 shows that this condition holds if λ 1 > −k 3 , which is why this condition is included in the second row. For the third row, we observe that λ ∞ = −k 3 (1 + h 0 ) < −k 3 , so the range 3 < < −λ 1 is always empty if λ 1 < λ ∞ .
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 k 1 ≤ k 4 and k 1 > k 4 . For the third row, we have used that λ 1 (k 1 , k 2 , k 4 ) is monotonic decreasing in λ 1 . Hence if k 1 ≤ k 4 , then λ 1 (k 1 , k 2 , k 4 ) ≥ λ 1 (k 4 , k 2 , k 4 ) = −k 4 (see (60)). Thus if λ 1 < λ ∞ = −k 3 (1 + h 0 ) and k 1 ≤ k 4 , then k 3 < − λ 1 1+h 0 ≤ k 4 1+h 0 ≤ k 4 1+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 λ ∞ = −k 3 (1 + h 0 ) and λ + * = −k 3 f + * (h 0 ).
We now consider a number of special cases that fill in the gaps that allow us to go from Tables 2 to 3.
• λ − * < λ 1 < 0, = 2 . When = 2 , then λ 1 = λ 4 and so there are two repeated eigenvalues that are the least negative. With some calculations, it can be seen that the eigenvectors v 1 and v 4 collide (the denominator (k 3 + λ 1 )( 3 + λ 1 ) + 2 k 3 h 0 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 v 4 = (0, 0, k 3 , k 3 + λ 4 ) and generically solutions will align with this eigenvector. In the proof of Lemma B.1, it is shown that k 3 +λ 4 > 0 when ≤ − 1 . Now 2 < − 1 for λ 1 > λ − * by (61) and this implies that the third and fourth components of v 4 have the same sign. The proof in Lemma B.1 now applies to prove that generically rebound occurs in this case.
• λ 1 = λ − * , = 2 = 3 = − 1 . In this case there is triple eigenvector with one eigenvector v 1 = v 4 = v 3 = (0, 0, k 3 , k 3 + λ 4 ) 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.
• λ + * < λ 1 < λ ∞ , = 2 with k 3 > k 4 1+m . This is similar to the first case as the two eigenvalues λ 1 and λ 4 coincide when = 2 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.
• λ + * < λ 1 < λ − * , = 3 with k 1 > k 4 or k 3 > −λ 1 . In this case, the real part of the complex pair λ 3,4 equals the eigenvalue λ 1 . Approaching the steady state via the complex pair leads to rebound. The proof of Lemma 5.14 also applies with = 3 provided that 3 + λ 1 < 0 and this occurs if k 3 > −λ 1 and so rebound occurs if the steady state is approached tangent to v 1 . 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.
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 λ 1 = λ + * intersects the line k 3 = k 4 1+m . Since lim k 1 →∞ = −(k 2 + k 4 ) (see (60)), if k 2 < k 4 f + * −(1+m) 1+m then the curve λ 1 = λ + * is always below k 3 = k 4 1+m . Also, as λ 1 (k 4 , k 2 , k 4 ) = −k 4 (see (60)), then it can be shown that the line k 3 = k 4 1+m is above the line λ + * at this point, and so the line can never go below the curve as the curve is monotonically increasing with k 1 .
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 0 > 0 such that rebound occurs for all satisfying 0 < < 0 (see Fig. 9). Thus, for all parameter values, rebound occurs for sufficiently small > 0. This might seem to contradict the result for the limit = 0 (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 x that represents the ligand and product in the compartments. The generalised TMDD model with feedback becomesẋ where x ∈ R n , n ≥ 1 and y, w ∈ R. We note that if f 1 depends only on x but not on y, or if f 2 depends only on y but not on x, then the coupling between x 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 f 1 and f 2 : • There exists an x 0 ∈ R n such that the Eqs. (63)-(65) have the steady state solution x = x 0 , y = w = 1. This implies that f 1 (x 0 , 1) = 0 and f 2 (x 0 , 1) = 0. 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., f 1 (x 0 , y) = 0 for all y ≥ 0 which implies that ∂f 1 ∂ y (x 0 , 1) = 0. Similarly, the function f 2 is not affected by the receptor dynamics if the ligand and product are at baseline, i.e., f 2 (x 0 , y) = 0 for all y ≥ 0 which implies that ∂ f 2 ∂ y (x 0 , 1) = 0. In fact, in both cases, all that we really require is the weaker derivative condition.
• The initial conditions are given by for some x 1 = x 0 with x 1 such that f 2 (x 1 , 1) < 0. Hence initially the amount of receptor will decrease. • The Jacobian matrix D x f 1 (x 0 , 1) has no zero or purely imaginary eigenvalues.
With these assumptions, the Jacobian matrix evaluated at the steady state is given by Note that the condition ∂f 1 ∂ y (x 0 , 1) = 0 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 D x f 1 (x 0 , 1) and the eigenvalues of the 2 × 2 bottom right matrix in J (x 0 , 1, 1). The global asymptotic stability of the steady state (x 0 , 1, 1) implies that all eigenvalues have negative or zero real part. Note that the eigenvalues of D x f 1 (x 0 , 1) do not depend on and that our final assumption implies that the eigenvalue closest to zero will have strictly negative real part. The 2 × 2 bottom right matrix is the same as in the basic TMDD model with feedback. We denote its eigenvalues by λ n and λ n−1 , where the eigenvalue λ n corresponds to λ 4 earlier and λ n−1 to λ 3 . Thus λ n is the closest eigenvalue to zero for sufficiently small .
We define˜ where˜ 2 is the value of at which the eigenvalue λ n coincides with the eigenvalue closest to zero of the matrix D x f 1 (x 0 , 1). Recall that − 1 is the value of for which λ n and λ n−1 become a complex pair of eigenvalues. We then have the following result. Theorem 6.1 If the above assumptions on f 1 and f 2 hold, then rebound occurs generically for 0 < <˜ 0 .
Proof The condition f 2 (x 1 , 1) < 0 ensures thatẏ(0) < 0, 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 2 × 2 matrix in J (x 0 , 1, 1) 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 by Non-dimensionalising as before, we obtain a new equation for the feedback variable w (= F) given byẇ = g(y, w), where g(y, w) = G(R 0 y, w) and 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 y = w = 1. 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 w = 1 of (67), assuming that y = 1 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 Proof We note that assumptions (h1 ) and (h2 ) imply that g(y, 1) > 0 for y < 1 sufficiently close to 1. The extra assumption requires that this holds for all y ∈ [0, 1). 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,ẇ(0) = 0 as before and in this case, we havë and so Lemma 5.9 holds. For Lemma 5.10, we note thaṫ w| w=1 = g(y, 1) > 0 using our stated assumption. Combined with Lemma 5.9, this implies that w(τ ) cannot touch or cross the line w = 1 when y(τ ) < 1. The second part of the proof of Lemma 5.10 shows by contradiction that y(τ ) = w(τ ) = 1 is impossible and the same proof holds in this case with h 0 replaced by g 1 . The Jacobian matrix associated with Eqs. (49) and (67) (analogous to (51)) is given by which has eigenvalues λ 3 = −k 3 , λ 4 = 0 when = 0 as previously. These two eigenvalues collide when the discriminant of the characteristic polynomial is zero, and this first occurs when With this definition of − 1 , we claim that Lemma B.1 holds also. To see this, we note that the eigenvector v 4 given by (59) is the same in this case, but with λ 4 as the eigenvalue of (69) that is closest to zero. We now have It is easily verified that k 3 > − 1 g 2 and this implies that k 3 + λ 4 > 0 for 0 < < − 1 . Thus, if 0 < < − 1 , then the third and fourth components of the eigenvector v 4 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 v 3 in (59), we note that it is easily verified that k 3 + λ 3 > 0 and so the same arguments apply in this case. The argument for the eigenvector v 2 is unchanged, while for v 1 , the requirement for the third and fourth components of the eigenvector to be the same is g 2 + λ 1 < 0 and it is easily verified that this holds. If = 2 corresponds to the point at which the eigenvalues λ 1 = −k 4 and λ 2 coincide, then the conditions of Lemma B.1 hold with 0 = min( − 1 , 2 ) and so we conclude that rebound therefore occurs for all satisfying 0 < < 0 .
This result shows that for this generalised model, rebound always occurs for sufficiently small > 0 even if there is no rebound without feedback ( = 0), as is the case for the original TMDD model with feedback moderator.
Clearly, if the bottom right 2 × 2 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 k e(L) = 0.024 day −1 , k e(P) = 0.201 day −1 , k out = 0.823 day −1 , R 0 = 2.688 nM, k off = 0.900 day −1 , k on = 0.592 (nM day) −1 , k in = k out R 0 = 2.212224 nMday −1 , and initial drug dose L 0 = 14.8148 nM. These parameter values have k e(L) < k e(P) 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 k e(L) < k e(P) 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 α > 0, even though there is no rebound at α = 0. The parameters for this example give λ 1 = −0.084, λ 4 = 0.126, k 3 = 0.517.  Fig. 9 and that rebound will happen for α < 0.135 (using α = /0.625). 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 H 0 = 1. The magnitude of the rebound is sufficiently small in this case that a nonlinear function H 1 (R) 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 α > 0 and sufficiently small. For this example, the rebound ends at approximately α = 0.98. In this case, a local maximum for R persists for larger values of α, but with R max /R 0 < 1. The value of time t max at which the maximum rebound occurs is also shown in Fig. 10. We note that as α → 0 then t max → ∞ 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 α = 1 since a local maximum persists in this case but at values lower than baseline, as already discussed.
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) are where X sc , X 1 and X 2 are the amounts of efalizumab in the depot, central and peripheral compartments respectively, X 3 is the total %CD11a on the surface of each T cell and X 4 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 X 1 X 3 K mc V c +X 1 instead of X 1 K mc V c +X 1 in the PK equations and have a more TMDDlike 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, X 3 is the total receptor, free and bound, and X 4 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 X = [X sc , X 1 , X 2 ], Y = X 3 and F = X 4 /k 03max . The differential equations for Y and F are then given by where The 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 by where F 0 is the positive root of the quadratic equation Substituting F = 1 + δ F into (72) gives a quadratic equation for δ F which has two negative roots, and this implies that F 0 ∈ (0, 1). It is then easily verified that all the conditions stated in Sect. 6.1 on the functions F 1 and F 2 hold. The particular function H that is used in this model is given by We non-dimensionalise, considering only Y and F, by defining the new variables The non-dimensional form of equations (70)-(71) is theṅ where dot denotes differentiation with respect to τ , x is the (unspecified) nondimensional form of X, and Note that we have used (72) in deriving h(y). All of our assumptions (h1)-(h4) on the function h(y) are now satisfied and 1). 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 k off from the stated value of 0.00154 day −1 to 0.0154 day −1 (the smaller value of k off 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 (0.033/150,000) × 10 −3 M = 2 × 10 −10 M, or 200pM. If we use the reported k off of 0.00154 day −1 and assume a typical k on (see Aston et al. 2011) of 0.592 nM −1 /day, the ratio k off /k on ≈ 2pM, i.e. 100-fold more potent than reported. A ten-times higher k off 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 k a = 0.242 day −1 , k 10 = 0.114 day −1 , k 12 = 0.097 day −1 , k 21 = 0.193 day −1 , k 30 = 0.444 day −1 , V c = 64.3 mL/kg, V m = 26.9 µg/mL, V m2 = 2.16 day −1 , K mc = 0.033 µg/mL, F a = 0.564, k 03max = 334 %CD11a/day. The plot of the total receptor X 3 in this case is the solid blue line shown on the left in Fig. 11. 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 X 4 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.

Efalizumab, no feedback
Total %CD11a Free %CD11a 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 Next we will interpret these results using our analysis. Using the parameter values in Ng et al. (2005), but the modified value of k off , we find that = 7.13 × 10 −3 , and 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 2×2 diagonal block. This is indeed the case since the eigenvalue closest to zero of this block is −1.37 × 10 −2 while the eigenvalue closest to zero from the remaining 3 × 3 block D x f 1 (0, 1) is −8.87 × 10 −2 . Thus, in the notation of Theorem 6.1, we have 0 < < 0 and so this Theorem says that we generically have rebound in the total receptor X 3 , as confirmed in Fig. 11 and by the numerical results in Ng et al. (2005).
The two eigenvalues of the lower 2 × 2 diagonal block collide at λ = −0.122 when = 0.0390, but since there are eigenvalues of D x f 1 (0, 1) 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 (k off ), as described in Lemma B.3.
Finally, we recall that the variable X 3 in the above model is the total of the free and bound receptor. The free receptor on T cells is given by K mc V c X 3 /(K mc V c + X 1 ) (Ng et al. 2005) and we note that if X 1 = 0, then the free receptor is the same as the total receptor. Now X 1 is the amount of efalizumab in the central compartment and while this is present, the total receptor X 3 is kept at a low level. However, once X 1 has been almost depleted, then X 3 starts to rise again and rebound occurs. Thus, X 1 is very low in the phase of the dynamics in which rebound in X 3 occurs, and this implies that the free receptor is very similar to X 3 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 (k e(P) < k e(L) and k e(P) < k out )).
At first sight, it might be surprising that there is always rebound for a very slow feedback response rate (0 < α 1), while there are large regions without rebound when the feedback moderator does not change and there is no feedback (α = 0). However, it can be seen that, in the absence of feedback dynamics (α = = 0), a larger initial condition for the moderator w will lead to a steady state with y = w 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 (α, > 0), 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 w > 1, 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 (α 1), 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 k e(L) , k out and k e(P) . An extra feature now is that the characteristics of the feedback function H (R) play a role as well if k e(L) > k e(P) . 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 (k e(L) < k e(P) ). 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 (k e(P) < k e (L) and k e(P) < k out ). However, if k out < k e(P) < k e(L) , 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 k on and k off 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 k on 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 k e(L) < k e(P) or small k out 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.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Proof of Lemma 4.2 Setting the derivatives to zero in (16) and (17) and solving for x and z, the only valid solution possible is x = z = 0. Substituting these into (18) with the derivative set to zero, we obtain By assumptions (h2) and (h3), y = 1 is the unique non-negative solution of this equation. Thus the unique steady state in the positive octant of Eqs. (16)-(18) is given by x = z = 0 and y = 1.
To show the global stability, we note that u = x +μz again decreases monotonically to zero using the same approach as in Lemma 2.3 of Aston et al. (2014) and since x and z must be non-negative, we conclude that lim τ →∞ x(τ ) = lim τ →∞ z(τ ) = 0, as for the original TMDD model. Similar to the original TMDD model, the y-axis, i.e, the manifold x = 0, z = 0, forms an invariant manifold asẋ | x=0, z=0 = 0 anḋ z | x=0, z=0 = 0. On this invariant manifold, the y-dynamics is given bẏ From the assumptions on the function h, it follows immediately that for 0 < y < 1, h(y) − y > 1 − y > 0; for y > 1, h(y) − y < 1 − y < 0; and for y = 1, we have h(y) − y = 1 − 1 = 0. Thus y = 1 is an attractor for the positive y-axis. Intuitively, these two observations together give that the steady state (x, y, z) = (0, 1, 0) is a global attractor for the positive octant. To make this formal, we consider the y-dynamics in a neighbourhood of the y-axis.
Thus for τ > T , the y-dynamics (18) can be estimated by If y(T ) > 1, then y(τ ) ≥ 1 on some τ -interval beyond T . In this interval we have h(y(τ )) ≤ 1, thuṡ Similarly, if y(T ) < 1, then y(τ ) ≤ 1 on some τ -interval beyond T . In this interval we have h(y(τ )) ≥ 1, thuṡ If y(τ ) > 1 for all τ ≥ T , then the estimate (73) is valid for all τ > T and we can conclude that there is some T 1 > 0 such that 1 ≤ y(τ ) ≤ 1 + 2 for τ > T 1 . Similarly, if y(τ ) < 1 for all τ ≥ T , then the estimate (74) is valid for all τ > T and we can again conclude that there is some T 1 > 0 such that 1 ≥ y(τ ) ≥ 1 − 2 for τ > T 1 However, it is also possible that there is some T ≥ T such that y( T ) = 1. In this case, (73) gives that if y(τ ) ≥ 1 in some interval [ T , T 1 ], then we can estimate for τ in this interval (74) gives that for τ in this interval, we have a, b) such that y(a) = 1 = y(b) and y(τ ) − 1 is either negative or positive for τ ∈ (a, b). On each of those intervals, one of the two estimates above can be used and hence we get that |y(τ ) − 1| < < 2 for all τ ≥ T . Combining these three cases, we can conclude that for all > 0, there is some T > 0 such that |y(τ ) − 1| < < 2 for all τ > T . Thus lim τ →∞ y(τ ) = 0 and we can conclude that all trajectories must converge to the steady state (x, y, z) = (0, 1, 0) as t → ∞.
Proof of Lemma 5.2 From (53) we note that Taking the limit in this expression as → ∞ gives the stated result.
Proof of Lemma 5. 3 We define the Lyapunov function For the steady state to be globally asymptotically stable, we have to prove that V is positive definite and continuous and thatV is negative definite and continuous relative to the steady state y = w = 1 on the whole of the phase space (Jordan and Smith 2007). By hypothesis (h4), V andV are continuous. To show that V is positive definite, we note that V (1, 1) = 0 and that hypothesis (h3) implies that We note that (75) implies that (1 − h(y))(y − 1) > 0 if y = 1 and soV < 0 if (y, w) = (1, 1). Also,V (1, 1) = 0. Thus,V is negative definite. As all the above conditions are satisfied, we conclude that the steady state (y, w) = (1, 1) is globally asymptotically stable.
Proof of Lemma 5. 4 We first observe that the nullclines for Eqs. (49), (50) are given by the two curves w = y, w = h(y) which intersect at the steady state y = w = 1. These two curves divide the phase plane into four regions in which the signs ofẏ andẇ do not change, as illustrated in Fig. 12.
We note from (55) that and it follows from this that We next observe from (52) and (53) that and it is easily verified from this that We now consider the two eigenvectors associated with the eigenvalues λ 3 and λ 4 , which are given by (54). The slopes of the linearised manifolds E 3,4 associated with the eigenvectors v 3,4 in the (y, w) plane are (k 3 + λ i )/k 3 , i = 3, 4. We now consider separately the two ranges in for which the eigenvalues are real.
If 0 < < − 1 , then from (76), we also have that < k 3 . From (77), this implies that λ i + k 3 > 0, i = 3, 4 and so the slopes of E 3 and E 4 are both positive and as λ 3,4 < 0, the slopes are also both less than one. Since λ 3 < λ 4 , the manifold associated with λ 4 has the steepest slope. This situation is shown in Fig. 13a.
Since λ 4 is the eigenvalue closest to zero, a generic trajectory will approach the steady state tangent to the linearised manifold E 4 . Now the global one-dimensional stable manifolds which are tangent to E 3 and E 4 at the steady state are invariant, and hence no trajectory can cross them. Thus, all initial conditions with y(0) < 1 but above the global stable manifold tangent to E 3 will rotate round the fixed point in a clockwise direction, as indicated in Fig. 12, and will approach the steady state tangent to E 4 from the right of the steady state, for which y > 1. Now this global manifold must itself spiral in to the steady state in a clockwise direction as shown in Fig. 12 and so cannot cross the line w = 1, 0 ≤ y ≤ 1, at least not unless it circles the steady state first. Thus, all initial conditions satisfying w(0) > 1, y(0) < 1 will give rise to trajectories with rebound as they approach the steady state, as will initial conditions with w(0) < 1, y(0) < 1 that lie above the global manifold that is tangent to E 3 at the steady state. When > + 1 , an analogous process can be followed to show that E 3 and E 4 both have negative slopes which are less than h (1) and this is shown in Fig. 13b. Again, any initial conditions to the right of the global manifold tangent to E 3 will rotate around the steady state in a clockwise direction and will approach the steady state with y > 1, giving rebound. As → ∞, we have already seen that λ 3 → −∞, and in this case, E 3 becomes vertical and so the region of initial conditions that gives rise to rebound, near to the steady state at least, becomes smaller as increases. Thus, for all ∈ ( + 1 , ∞), there are always some initial conditions with w(0) > 1, y(0) < 1 which will give rise to rebound, as claimed.
We note that as the global stable manifolds may spiral in to the steady state in a clockwise direction, we cannot say whether initial conditions satisfying w(0) < 1, y(0) < 1 are above or below this global manifold, and so we cannot comment on the existence or non-existence of rebound in this case. When = ± 1 , the eigenvalues λ 3 and λ 4 collide as well as their eigenvectors. Instead of two eigenvectors, there is now one eigenvector v 3 = v 4 and a generalised eigenvector. All solutions will approach the steady state via the eigenvector, see Jordan and Smith (2007). The analysis for the approach to the steady state for the case 0 < < − 1 can be extended to the case = − 1 and similarly, the case > + 1 can be extended to = + 1 . The last case to consider are the complex eigenvalues when − 1 < < + 1 . When is in this range, the eigenvalues λ 3,4 are a complex conjugate pair. Thus close to the steady state, trajectories will spiral around it. Each time the trajectory loops round the steady state, there will be a part of the trajectory with y > 1 and so rebound occurs infinitely often, but with exponentially decreasing magnitude, which is determined by Re(λ 3,4 ) = −(k 3 + )/2.
Proof of Lemma 5.5 The phase space is bounded by the planes x = 0, y = 0, z = 0 and w = 0. As in the proof of Lemma 4.1, we have to show that the derivative normal to each of these planes is non-negative (Smith 1995).
The proof that the derivative normal to the planes x = 0 and z = 0 is non-negative is the same as in the proof of Lemma 4.1. On the plane y = 0, we find from (22) thaṫ y| y=0 = k 3 w + k 2 z and soẏ| y=0 ≥ 0 when w, z ≥ 0. Similarly on the plane w = 0, from (23), we obtaiṅ w| w=0 = h(y) and soẇ| w=0 > 0 for all y ≥ 0 using assumption (h1).
Proof of Lemma 5.6 Solving the equationsẋ =ẏ =ż = 0 for x, y and z again gives two solutions, but one of these has y negative, so is not physically relevant and is not in the positive part of the phase plane considered in this lemma. The other solution is given by Substituting for y into the equationẇ = 0 gives As in the proof of Lemma 4.2, we see that w = 1 is a solution of this equation using (h2) and assumption (h3) ensures that this solution is unique.
Proof of Lemma 5.7 We showed in Aston et al. (2014) that λ 1 and λ 2 are always real and negative. Furthermore, Lemma 5.2 gives that λ 3 and λ 4 are always negative and real or complex with negative real part. So the steady state is linearly stable.

Appendix B Proof of rebound lemmas
Lemma 5.13 is proved by assuming that λ 4 is closest to zero and then considering the three intervals for which the value of λ 4 is qualitatively different.
Lemma B.1 If 0 < ≤ − 1 and the eigenvalue λ 4 is closest to zero then rebound occurs.
Proof First we consider the case 0 < < − 1 . Generically the trajectory must approach the steady state tangent to the eigenvector v 4 given in (59) since λ 4 is real and is the eigenvalue closest to zero. We note that We see from (53) that the term under the square root is positive if λ 4 is real (i.e., / ∈ ( − 1 , + 1 )), and thus the expression above immediately implies that k 3 + λ 4 > 0 if < − 1 < k 3 and k 3 +λ 4 < 0 if > + 1 > k 3 . Therefore, if 0 < < − 1 , then the third and fourth components of the eigenvector v 4 have the same sign. We therefore have two cases to consider. In the first case, y and w both approach the steady state from above, and clearly rebound occurs in this case. Alternatively, y and w both approach the steady state from below and Corollary 5.11 implies that there is rebound in this case also.
As in the proof of Theorem 3.5 in Aston et al. (2014), we must again consider the possibility that for particular parameter values, the generic situation considered above does not occur, and so the trajectory approaches the steady state tangent to one of the other eigenvectors.
We first consider the eigenvector v 3 (see (59)). For this vector, we note that k 3 + λ 3 > 0, by a similar argument to the one above for k 3 + λ 4 , and so the third and fourth components of the eigenvector v 3 have the same sign. Thus, the above proof by contradiction also holds if the trajectory approaches the steady state tangent to v 3 .
For the eigenvector v 2 , we note, as in the proof of Theorem 3.5 of Aston et al. (2014), that this eigenvector points out of the phase space, since the first and second components have opposite sign, and so no trajectory can approach the steady state tangent to this eigenvector.
Finally we consider the eigenvector v 1 . We note that the third and fourth components of this eigenvector will have the same sign if + λ 1 < 0. Since λ 4 is real and closest to zero, the definition of λ 4 (53) gives + λ 1 < + λ 4 = 1 2 − k 3 + (k 3 − ) 2 − 4 k 3 h 0 < 0 as < − 1 < k 3 . So we can conclude that the third and fourth entries in the eigenvector v 1 have the same sign, and so the above proof by contradiction also holds if the trajectory approaches the steady state tangent to v 1 .
This covers all possible cases for < − 1 , and so we conclude that rebound must occur for those values.
Finally we consider the case = − 1 , which is when the Jacobian matrix has two repeated real eigenvalues λ 3 = λ 4 (= λ − * ) that are the least negative. In this case, the two corresponding eigenvectors satisfy v 3 = v 4 and there is a generalised eigenvector v g which is given by v g = (0, 0, 0, 1) The solutions on the two-dimensional linearised manifold spanned by v 3 and v g are linear combinations of e λ 3 t (v g + tv 3 ) and e λ 3 t v 3 (Jordan and Smith 2007) and asymptotically will align with the eigenvector v 3 . It is easily verified that k 3 + λ − * > 0, and so the third and fourth components of v 3 have the same sign. The same arguments as above show that rebound must then occur in this case, and that the trajectory cannot approach the steady state tangent to v 1 or v 2 .
Next we consider the case that > + 1 (k 3 ) and λ 4 is real and the eigenvalue closest to zero. This can only happen if λ 1 < λ ∞ . In the case of the direct feedback, we analysed this region by using that the total amount of receptor (v = y + z) stays on one side of the plane v = 1, see Lemma 4.7. As observed earlier, this property does not hold any more for the feedback with a moderator, but we can define a related quantity and show a partial result.
At baseline, we haveṽ = 1+ k 3 and the initial condition is also in the planeṽ = 1+ k 3 . By using this plane, we will now prove that there is rebound if k 3 is sufficiently large compared to k 4 .
Since λ 4 is the eigenvalue closest to zero, generically the solution will approach the steady state along the line spanned by the eigenvector v 4 . Furthermore, the outward normal to the planeṽ = 1 + k 3 is ∇ṽ = (0, 1, 1, k 3 / ). Thus with the expression for v 4 in (59) and > + 1 , it follows that its fourth component is As the solution has y < 1 and w > 1, the solution comes in via −v 4 . But this vector is below the planeṽ = 1 + k 3 as as, by the definition of λ 4 , we have + k 3 + λ 4 = 1 2 (k 3 + + √ D) > 0. This contradicts the fact that the solution is always above the plane, hence our assumption that there is no rebound is false.
Finally, we consider the case = + 1 . This case is similar to that when = − 1 in the proof of Lemma B.1 in that there is a repeated eigenvalue λ 3 = λ 4 with one eigenvector v 3 = v 4 and a generalised eigenvector v g . The solution will again asymptotically align with the eigenvector v 4 which has fourth component k 3 + λ + * < 0 and so the above proof by contradiction works in this case also.
Next we look at the case that the eigenvalue λ 4 is complex (i.e, − 1 < < + 1 ) and its real part is closest to zero.

Lemma B.3
If the eigenvalue λ 4 is complex and its real part is larger than λ 1 , hence closer to zero, then y generically converges to its steady state value y = 1 in an oscillatory fashion, and so rebound occurs infinitely many times, but with exponentially decreasing magnitude.
Proof If the complex conjugate pair of eigenvalues λ 3,4 have real part that is closer to zero than λ 1 , then the Hartman-Grobman Theorem (Crawford 1991) gives that generically trajectories will converge to the steady state in the phase space tangent to the (y, w) plane and will spiral around the steady state (y, w) = (1, 1). Each time the trajectory loops round the steady state, there will be a part of the trajectory with y > 1 and so rebound occurs infinitely often, but with exponentially decreasing magnitude, which is determined by Re(λ 3,4 ) = −(k 3 + )/2.
As in the previous case, we also need to consider the possibility that for particular parameter values, the trajectory might approach tangent to a different eigenvector. We again note that the eigenvector v 2 points out of the phase space, and so a trajectory cannot approach the steady state tangent to this eigenvector. However, the eigenvector v 1 points into the phase space, since the first and second components of the eigenvector are both positive. Thus, we are unable to exclude the possibility that an orbit could approach the steady state tangent to this eigenvector for particular parameter values. However, this is very unlikely to occur in practice and is not a generic case.
Proof of Lemma 5.13 This follows immediately by combining Lemmas B.1-B.3 and using that − 1 < < + 1 implies that λ 4 is complex.
Next we will consider the case that λ 1 is the eigenvalue closest to zero. The eigenvector v 1 corresponding to the eigenvalue λ 1 is given in (59). The first two entries in this eigenvector are the same as we had previously for the case with no feedback, and since both of these quantities are positive (see Lemma 2.7 of Aston et al. 2014), we must take a positive multiple of v 1 to ensure that x and z are positive in the direction of this eigenvector. We now consider the signs of the third and fourth components of this vector, although we note that since the steady state values for y and w are non-zero, there is no requirement that these entries be positive. First we determine the sign of the denominator in the expressions for the third and fourth components of v 1 .
Proof of Lemma 5.14 The ratio of the third and fourth components of v 1 (which we denote by v 1,3 and v 1,4 respectively) is given by Thus, if + λ 1 < 0, then v 1,3 and v 1,4 have the same sign and the proof of rebound given in Lemma B.1 holds in this case also.
If k 1 > k 4 , then we note from Lemma 2.8 of Aston et al. (2014) that k 4 + λ 1 < 0 and so, using the result of Lemma B.4, we then see that v 1,4 < 0 in this case. Since we are taking a positive multiple of the eigenvector, this implies that w must approach its steady state value from below and so, by Corollary 5.11, rebound must occur in this case also.
As in previous cases, we should again consider the possibility that a trajectory could approach the steady state tangent to a different eigenvector. As previously, the eigenvector v 2 again points out of the invariant region since the first two entries have opposite sign. However, it is not possible to determine the global invariant manifold that approaches the steady state tangent to the eigenvectors v 3 or v 4 and so we can only state a generic result.