Laurent Reynaud
email: laurent.reynaud@orange-ftgroup.com

Tinku Rasheed
email: tinku.rasheed@create-net.org

Sithamparanathan Kandeepan
email: kandeepan@ieee.org

An Integrated Aerial Telecommunications Network that Supports Emergency Traffic

Keywords: aerial networks, emergency communications, low altitude platforms, high altitude platforms, controlled mobility, disaster medical assistance team, Hurricane Katrina I

This paper outlines how an aerial telecommunications network can optimally meet the stringent needs of emergency relief and recovery operations. We propose a novel architecture, made of an integrated and highly dynamic multi-purpose aerial telecommunications infrastructure that can be contextually extended with fast-deploying high or low altitude platforms. In particular, we analyze the interest and challenges of adapting core concepts from substitution networks and controlled mobility mechanisms, so that a base network can be seamlessly augmented, both in terms of capacity and functions. We give an estimation of the emergency traffic supported by the lower altitude platforms in an example scenario and discuss the challenges posed by this architecture, notably in terms of disaster resilience and ability to efficiently provide sustained first responder communications.

INTRODUCTION

In the aftermath of a disaster, whether natural or man-made, it is of the utmost importance that rescue teams and other potentially life-saving response operations can rely on an efficient emergency communications system. Unfortunately, terrestrial networks may be partially or totally compromised, at least during the early stages of the post-disaster period. Similarly, alternative communications systems may prove less than effective. As an illustration, the report of the Select Bipartisan Committee to investigate the response to Hurricane Katrina [START_REF] Davis | A Failure of Initiative: Final Report of the Select Bipartisan Committee to Investigate the Preparation for and Response to Hurricane Katrina[END_REF] clearly concluded that massive network failures and severe interoperability issues with the remaining communications systems greatly hindered rescue effort, situational awareness and command and control operations. For the support of reliable emergency communications, recovery solutions relying on satellite communications have long been investigated through multiple research activities including the European projects Real-Time Emergency Management via Satellite (REMSAT) [START_REF] Pace | Disaster monitoring and mitigation using aerospace technologies and integrated telecommunication networks[END_REF] and Emergency Management by Satellite Communications (EMERGSAT) [START_REF] Pace | Disaster monitoring and mitigation using aerospace technologies and integrated telecommunication networks[END_REF].

In this context, aircraft-based telecommunications, which we investigate through this paper, have also gained momentum in the recent years, with many types of aerial platforms seen as a possible extension or alternative to traditional equipments from terrestrial and satellite wireless networks. Those platforms, either aerostats or aerodynes, evolve at various altitudes and are named Low Altitude Platforms (LAPs) [START_REF] Hariyanto | Emergency broadband access network using low altitude platform[END_REF] or High Altitude Platforms (HAPs) [START_REF] Aragón-Zavala | High-Altitude Platforms for Wireless Communications[END_REF], [START_REF] Grace | Broadband Communications via High-Altitude Platforms[END_REF], whether they fly in the lower troposphere or in the stratosphere. Depending on multiple factors including their altitude, available power for the payload, the type of antenna and radio technology, an aerial platform can cover areas of up to 300 km radius. Moreover, many aerial vehicles are able to keep a quasi-stationary position. Direct Line Of Sight (LOS) conditions with terrestrial equipments are thus often met. As a consequence, many use cases that take benefit of aerial platforms have been recently envisioned [START_REF] Aragón-Zavala | High-Altitude Platforms for Wireless Communications[END_REF], [START_REF] Grace | Broadband Communications via High-Altitude Platforms[END_REF], including emergency communications. Due to their payload flexibility, relatively fast deployment times and low altitude (hence allowing reduced radio propagation delays), HAPs and LAPs are increasingly seen as a viable solution to support relief services. They can quickly restore communications in the affected zones without mandatory support of the terrestrial infrastructure, which may be at least partially inoperable. Moreover, the aerial platforms can be effectively deployed together with the existing satelliteterrestrial recovery infrastructure, to provide speedy deployment able to cover larger areas and to maintain effective communication links supporting adequate capacities and bandwidth demands of the terrestrial first responders. However, it is observed that relief scenarios based on aerial platforms lack a global architecture that can federate and support all the requirements of the multiple actors involved in emergency relief and recovery operations. On the opposite, each service or resource requirement related with an emergency role (e.g. search and rescue personnel also known as first responders, command and control operations and critical infrastructure repair teams) is currently addressed independently, in a standalone fashion.

In this paper, we present an integrated Aerial Telecommunications Network (ATN) architecture for emergency communications, constituting an integrated and highly dynamic multi-purpose aerial telecommunications infrastructure that can be contextually extended with fastdeploying aerial platforms. We select a set of typical emergency services needed by rescuers, and outline the subsequent impacts in term of required capacity on the considered network. Bandwidth requirements are estimated for each type of traffic. We then argue that such an integrated emergency communications platform is inevitable given the predictability of different disasters and the flexibility provided with fast and augmented critical infrastructure restoration support that can meet the variability of the estimated traffic requirements. We also analyze the interest and challenges of adapting core concepts from substitution networks and controlled mobility mechanisms, so that a base network can be seamlessly augmented in terms of coverage, capacity and functions.

II. RELATED WORK

Several studies identified the benefits of aerial platforms with the objective of supporting wireless communications in emergency situations. In these works, targeted airborne platforms are not always designed for higher elevations. LAPs in particular are not only seen as a fast and convenient way to experiment telecommunications payloads with relatively inexpensive aerial vehicles, but can also address actual scenario requirements where limited coverage, due to low altitude, is acceptable. The authors in [START_REF] Hariyanto | Emergency broadband access network using low altitude platform[END_REF] and [START_REF] Qiantori | An Emergency Medical Communications System by Low Altitude Platform at the Early Stages of a Natural Disaster in Indonesia[END_REF] sought to evaluate several LAP vehicles for relief support in several regions of Indonesia particularly exposed to a large number of potential threats, including droughts, floods, landslides, earthquakes and volcanic eruptions. In this context, a simple tethered balloon, flying at a maximum altitude of 440 m, can carry an IEEE 802.11a/g payload and cover an area of about 72 km 2 . Naturally, in these cases, the covered area and maximum system throughput are relatively limited and only LOS communications are supported. Other radio technologies were also investigated, such as in [START_REF] Zhu | High altitude platform-based two-hop relaying emergency communications schemes[END_REF]. This work, outlining a multihop architecture using relaying features specified in the IEEE 802.16j standard, is not specific to HAPs and may also be considered for LAPs. Still, no performance evaluation was presented and the study thus remains speculative in this regard. In contrast, a practical approach was used in [START_REF] Huang | Emergency Communication for Electric Power System Based on Airship[END_REF] where emergency communications for electric power systems were investigated. It was shown that, compared to satellite systems, HAPs could favorably meet the stringent needs of power grid security in terms of both bandwidth and especially delay.

However, in a more general case, the support of emergency communications is investigated at higher altitudes. In [START_REF] Pace | Disaster monitoring and mitigation using aerospace technologies and integrated telecommunication networks[END_REF], the use of satellites to recover communications in disaster-affected areas is investigated, and several deployment scenarios that combine the features of low and geostationary earth orbiting satellites and HAPs are proposed. Although the multiple usecases clarify the advantages of integrated aerial infrastructures in disaster relief situations, no insight is given as to when and how the considered infrastructures are set up. The author of [START_REF] Deaton | High Altitude Platforms for Disaster Recovery: Capabilities, Strategies, and Techniques for Emergency Telecommunications[END_REF] investigated how first-responder and emergency answering services such as 911 or 112 could be supported by HAPs. One of the major considerations in this paper was to propose a monolithic HAP payload architecture that would enable interoperability with most representative first-responder technologies. Also, as cellular networks constitute an important communications asset for the various rescue and relief teams, the support of several technologies may be required. As a result, the specified payload consists of multiple access technologies with adapted antennas, Voice over IP (VoIP) and radio gateways, a full-featured IP core and several adapted services. It is thus unrealistic to consider that this payload architecture could be deployed according to its full specifications on a currently available HAP airplane or airship, without severe capacity limitations for each of the supported radio technologies.

III. A FLEXIBLE AERIAL TELECOMMUNICATIONS NETWORK

In our architecture, an ATN is a collection of aircrafts able to communicate with other nodes through adapted payloads. Compared to existing network topologies which include airborne platforms [START_REF] Aragón-Zavala | High-Altitude Platforms for Wireless Communications[END_REF], [START_REF] Grace | Broadband Communications via High-Altitude Platforms[END_REF], [START_REF] Pace | Disaster monitoring and mitigation using aerospace technologies and integrated telecommunication networks[END_REF], [START_REF] Deaton | High Altitude Platforms for Disaster Recovery: Capabilities, Strategies, and Techniques for Emergency Telecommunications[END_REF], we defined a multi-level hierarchical aerial topology. These levels, even if it is the case in this paper, are not necessarily related to different ranges of altitude. Instead, they are primarily intended to identify several sets of airborne platforms with common key properties such as mobility patterns, payload type, platform capacity and supported radio technologies. As an application of this general principle, Fig. 1 illustrates the desired network architecture for disaster relief operations, with two distinct levels:

• The higher level of aerial platforms consists of one or several quasi-stationary vehicles, integrated to satellite and terrestrial networks via high-capacity links. Naturally, the availability of terrestrial backhaul links also depends on the extent of damage caused by the considered disaster. However, remote access to ground stations which are located out of the damage zone should be facilitated by the support of high capacity Inter-Platform Links (IPL) and multi-hop relay mechanisms. These nodes may already be in place prior to the disaster, provided a commercial aerial network is already operated on the zone or the disaster can be anticipated to some extent. Moreover, they are assumed to be evolving at high altitudes, in order to benefit from HAPs key properties such as large terrestrial coverage, contained average wind intensity and favorable radio propagation conditions between both terrestrial and satellite segments.

• The lower level of aerial platforms is made of all the nodes that are set up and deployed for specific recovery missions. These platforms may have different mobility patterns, flight durations, payload types and communication links. As depicted by Fig. 1, the communications in this heterogeneous environment are facilitated by the support of IPL with neighboring same-level platforms, but also with higher level platforms. Note that IPL from different levels do not necessarily use the same radio technologies and the same types of antennas. In fact, higher level IPL are expected to use directive antennas and Free Space Optics (FSO), because of the relative stability of the vehicles and the favorable radio propagation conditions. However, due to the more dynamic mobility patterns and design constraints related to the lower level vehicles (including cost and weight), sectoral and omnidirectional antennas and other radio technologies are generally preferred for higher level IPL.

As a result, the ATN can be hierarchically partitioned into clusters [START_REF] Rasheed | A stable clustering scheme for large scale mobile ad hoc networks[END_REF], as illustrated by Fig. 2. This topological representation offers many advantages. In particular, an adapted large scale routing protocol could inherently differentiate the types of links (backhaul links and IPL) by using routing schemes adapted to each partitioned zone of the network. Furthermore, our ATN architecture is also designed to use the concept of controlled mobility [START_REF] Rao | Controlled Node Mobility: A Mechanism for Energy Load Distribution in Sensor Networks[END_REF], as Fig. 3 illustrates. In essence, a base network, experiencing troubles such as communication link instability or congestion, may be temporarily supplemented with additional wireless dirigible nodes. Moreover, an integrated ATN which federates multiple types of emergency recovery equipment is complex and costly to deploy, with a high demand in term of ground staff dedicated to the network operations [START_REF] Grace | Broadband Communications via High-Altitude Platforms[END_REF]. Controlled mobility is intended to alleviate this demand with self organization mechanisms that take benefit of the unique mobility capabilities of aerial platforms. This way, emergency personnel may be partially relieved of some remote topology control and vehicle guidance tasks and be allowed to concentrate on their respective core rescue and recovery missions. The envisioned way to enforce controlled mobility within the integrated ATN is manifold. Yet, two approaches in particular can be outlined:

• The base network may be made of LAPs (each node following a predetermined mobility pattern) and the terrestrial infrastructure equipment left undamaged by the disaster. In this scenario, controlled mobility is applied to the HAPs, which are moved so that they can provide maximum coverage and bandwidth to terrestrial and LAP equipments.

• This base network may instead consist of the HAPs, with quasi-stationary positions, plus, again, undamaged or repaired terrestrial equipment. Controlled mobility is here applied to the LAPs. LAP default mobility patterns will thus be altered so that they can extend the base network's capacities where it is best needed. In the rest of this paper, we consider and discuss this approach.

One of the major expected benefits of the controlled mobility and other discussed ATN mechanisms is the ability to meet applicative demands in a given area with a relatively low number of nodes. In the next section, we seek to illustrate the corresponding ideal case, where a minimum number of LAPs can optimally support a variable emergency traffic demand.

IV. EMERGENCY TRAFFIC REQUIREMENTS

Out of the various services used by the emergency response teams, we chose to estimate the global voice emergency and medical care video traffics. In order to perform this estimation, we contextualized the emergency traffic with parameters consistent with the 2005 Hurricane Katrina disaster. Among those parameters are the affected area (233,000 km 2 [START_REF] Qiantori | An Emergency Medical Communications System by Low Altitude Platform at the Early Stages of a Natural Disaster in Indonesia[END_REF]), the extent and duration of the base mobile communication network outage [START_REF] Deaton | High Altitude Platforms for Disaster Recovery: Capabilities, Strategies, and Techniques for Emergency Telecommunications[END_REF] and the estimation of the population of the first responders and related traffic [START_REF] Deaton | High Altitude Platforms for Disaster Recovery: Capabilities, Strategies, and Techniques for Emergency Telecommunications[END_REF], [START_REF] Davis | A Failure of Initiative: Final Report of the Select Bipartisan Committee to Investigate the Preparation for and Response to Hurricane Katrina[END_REF]. Note that we took the example of the Hurricane Katrina disaster because this largescale event is well-documented, with in particular abundant literature about the composition and actions of the different response teams and about the coordination problems (partly due to inadequate communication systems) experienced by rescuers.

A. Voice emergency traffic estimation

We define the voice emergency traffic as the total number of emergency voice calls passed by the public and first responders (e.g. rescuers, army, emergency medical care, nonmedical care and public safety teams). We base our estimation on the related traffic pattern estimated in [START_REF] Deaton | High Altitude Platforms for Disaster Recovery: Capabilities, Strategies, and Techniques for Emergency Telecommunications[END_REF] for the 20 days following Hurricane Katrina landfall, and similarly perform a normalization of the daily traffic corresponding to a daylight rescue. For each day of the corresponding period, we calculate the busy hour traffic and take this conservative value as the traffic which the ATN must support. These values are plotted in Fig. 4. Furthermore, as our considered ATN architecture is intended to achieve voice communications interoperability through a VoIP system that supports adapted gateways, we set a series of plausible hypotheses about the used radio technologies and VoIP codecs and configuration. According to [START_REF] Qiantori | An Emergency Medical Communications System by Low Altitude Platform at the Early Stages of a Natural Disaster in Indonesia[END_REF], 78 VoIP calls can be simultaneously established with an acceptable quality from a LAP evolving at an altitude of 440 m and an achievable capacity of 54 Mb/s anywhere in a covered zone of 47.39 km 2 . In these conditions, with an 802.11g communication interface, an Adaptive Multi-Rate (AMR) codec at a rate of 12.2 kb/s and a sample period of 20 ms, a Mean Opinion Score (MOS) of 3.8 can be achieved. Hence, using these parameters for our estimation, we calculated and represented in Fig. 4 the minimum number of LAPs required to support the corresponding traffic. A maximum of 475 concurrent voice sessions supported by 7 LAPs is attained on day 15.

B. Video emergency traffic estimation

In our estimation, the video emergency traffic is emitted solely by the medical assistance teams, and more precisely, in the context of Hurricane Katrina, by the Disaster Medical Assistance Teams (DMATs). In the USA, DMATs, which belong to the National Disaster Medical Systems (NDMS), are teams of physicians, nurses, technicians, and support staff. Their main mission is to provide emergency care, perform triage and resuscitation, stabilize and prepare for evacuation. At the time of Hurricane Katrina landfall, there were 45 DMATs spread out across the USA. When deployed on a disaster zone, DMATs are typically made of 35 members [START_REF]Office of the Assistant Secretary for Preparedness and Response[END_REF].

According to [START_REF] Davis | A Failure of Initiative: Final Report of the Select Bipartisan Committee to Investigate the Preparation for and Response to Hurricane Katrina[END_REF], there were 9 DMATs initially prepositioned before Katrina landfall. Eventually, all the DMATs were progressively deployed in response to the disaster [START_REF] Lister | Hurricane Katrina: The Public Health and Medical Response[END_REF]. Our estimation reflects the ratio of prepositioned medical care units and the progressive activation of supplementary teams. Also, we took the following assumptions:

• Each deployed DMAT is made of 35 members, including 4 physicians. Those physicians work in pairs, and perform alternating 12-hour shifts. We thus consider that, at any time, there are two active physicians per DMAT. At any time, the number of video streams per DMAT is equal to the number of active physicians, i.e. 2.

• 9 DMATs are active at day 1, and this number of active DMATs progressively reaches 45 at day 15. Moreover, when deployed, we assume that each DMAT is fully operational for the remaining analyzed period.

• Also, we consider that with video sessions based on the H.264 codec with a bit rate of 384 kb/s, most medical video teleconsultation services can be supported with an acceptable quality [START_REF] Komnakos | Performance Evaluation of an Enhanced Uplink 3.5G System for Mobile Healthcare Applications[END_REF]. With the same configuration in terms of altitude, ground coverage and communication interfaces as used when estimating voice traffic, each LAP can support 18 simultaneous H.264 video sessions at the given bit rate.

With these hypotheses, we estimated and represented in Fig. 5 the number of active DMATs for the observed 20-day period. The subsequent number of required concurrent video sessions is easily deduced, and allows estimating the minimum number of LAPs (also represented in Fig. 5) able to support the corresponding traffic.

C. Aggregate emergency traffic estimation

While the required number of LAPs for aggregate voice and video traffic generally depends on whether the LAPs are dedicated to one or both traffic types, in our example however, the same maximum value of 12 LAPs is attained in both cases, at days 15, 16 and 18. Furthermore, it shall be noted that this estimation relies on the ability to serve every applicative demand with a relatively low LAP coverage. For instance, in our example, the maximum covered LAP area without overlap is 568.68 km 2 , and covers less than 0.25% of the overall disaster area. In the context of Hurricane Katrina, this low proportion is not unrealistic, considering the very low population density in a large part of the affected area, and the large number of rescue teams in a few key locations [START_REF] Davis | A Failure of Initiative: Final Report of the Select Bipartisan Committee to Investigate the Preparation for and Response to Hurricane Katrina[END_REF]. Nevertheless, some rescue teams may be mobile and the emergency traffic may vary with time and location. As a result, the sustained performance of the considered architecture strongly relies on the ability to take advantage of the HAP larger coverage and to enforce efficient ATN mechanisms, including controlled mobility. The different payload capabilities of the HAPs and LAPs will determine the number of modules that can be carried by these platforms and eventually the diversity of the services that can be provisioned from the ATNs. Moreover, platform station keeping, hand-off considerations even for fixed stations due to the platform movement and payload power bring lot of emphasis on the reliability of the aerial platforms in providing effective communications to the terrestrial stations.

The technological challenges and system issues that need to be addressed for effective integration of ATN platforms for emergency communications include the stratospherical environment and thermal conditions, launch and recovery, energy sources and platform aerodynamics. Apart from these, the technological choices for communication between the terrestrial responders and the aerial stations are an important aspect which governs the design of the system architecture, frequency planning, resource planning and network protocols. Moreover, the scale of disaster response requires flexible and incremental deployment strategies, including controlled mobility mechanisms and subsequently node coordination mechanisms with adequate positioning and clustering approaches to ensure maximum system availability and reliability.

Though energy efficiency is not the prime concern in a disaster scenario, it is of utmost importance to ensure the maximum operational lifetime for the ATN platforms and the terrestrial devices alike. This requires both efficient cooperative and context aware mechanisms which reduce the overhead in communications and techniques that scale down the power consumption in the uplink communications, apart from developing intelligent and energy efficient ground device interfaces. In addition, the ATN communications should complement well the established terrestrial disaster communication technologies like the TETRA wireless platforms and the COLT (Cells on Light Trucks) / COW (Cells on Wheels) platforms deployed by rescue operators [START_REF] Deaton | High Altitude Platforms for Disaster Recovery: Capabilities, Strategies, and Techniques for Emergency Telecommunications[END_REF]. The seamless integration of ATN with the satellite-terrestrial emergency service platforms and radical changes in the aeronautical and radio regulations can make integrated ATN a reality in providing emergency communications servicing large scale disasters.

VI. CONCLUSION

In this paper, we discussed a novel integrated aerial architecture for emergency communications. We selected a subset of typical voice and video emergency services needed by rescuers, and estimated the subsequent impacts in terms of bandwidth requirements through the example of Hurricane Katrina. Estimation results showed that, in a case where traffic is ideally shared on the aerial platforms, a small number of LAPs have enough capacity to support a variable emergency traffic emitted on a large area. However, their limited coverage being an issue in large scale events, we also argued that the sustained performance of this infrastructure, notably in term of coverage, strongly relies on the ability to take advantage of the multi-level topology, including higher altitude platforms, and to enforce dedicated mechanisms, including controlled mobility.

In the future, we intend to study the impact of emergency traffic on the various communication interfaces of the considered architecture. In particular, we are interested in determining what proportion of the emergency traffic is meant to remain local, regional or to be transmitted through the backhaul links. We will also seek to further evaluate the interest of controlled mobility to support efficient cooperative communications between the different aerial platforms.

Figure 1 .

 1 Figure 1. A 2-level aerial network adapted to disaster relief operations.Aerial vehicles are interconnected through IPL.

Figure 2 .

 2 Figure 2. Hierarchical representation of the multi-level ATN topology.

Figure 3 .

 3 Figure 3. The concept of controlled mobility. In case a, a communication link (represented by dotted lines) cannot meet the traffic demand. In case b, mobile nodes are introduced and placed so that an additional multi-hop route is established between the nodes with the mentioned link. If another link experiences troubles, as represented in case c, the mobile nodes can be moved accordingly.

  V. MAIN CHALLENGES FOR FUTURE INTEGRATED ATNOne of the prominent challenges affecting the development of integrated ATN's is the technological maturity of the platforms themselves[START_REF] Grace | Broadband Communications via High-Altitude Platforms[END_REF]. The airship-based HAPs and unmanned UAVs, both based on fuel and solar energy are at different levels of development. They are governed not only by technical and manufacturing constraints but also by both political and economic factors.

Figure 4 .

 4 Figure 4. Evolution of the numbers of emergency voice traffic and required LAPs to support the subsequent concurrent voice sessions (AMR codec at a rate of 12.2 kb/s with a sample period of 20 ms and a MOS of 3.8).

Figure 5 .

 5 Figure 5. Evolution of the numbers of active DMATs and required LAPs to support the subsequent H.264 video traffic.

This work was partially supported by the ANR RESCUE project, grant ANR-10-VERS-003 of the French Agence Nationale de la Recherche.


Philip J Aston
email: p.aston@surrey.ac.uk

Gianne Derks
email: g.derks@surrey.ac.uk

Adewale Raji

Balaji M Agoram
email: agoramb@medimmune.com

Piet H Van Der Graaf
email: piet.van.der.graaf@pfizer.com

Mathematical Analysis of the Pharmacokinetic-Pharmacodynamic (PKPD) Behaviour of Monoclonal Antibodies: Predicting in vivo Potency

Keywords: Pharmacology, efficacy, affinity, target-mediated drug disposition, IgE

We consider the relationship between the target affinity of a monoclonal antibody and its in vivo potency. The dynamics of the system is described mathematically by a target-mediated drug disposition model. As a measure of potency, we consider the minimum level of the free receptor following a single bolus injection of the ligand into the plasma compartment. From the differential equations, we derive two expressions for this minimum level in terms of the parameters of the problem, one of which is valid over the full range of values of the equilibrium dissociation constant K D and the other which is valid only for a large drug dose or for a small value of K D . Both of these formulae show that the potency achieved by increasing the association constant k on can be very different from the potency achieved by decreasing the dissociation constant k off . In particular, there is a saturation effect when decreasing k off where the increase in potency that can be achieved is limited, whereas there is no such effect when increasing k on . Thus, for certain monoclonal antibodies, an increase in potency may be better achieved by increasing k on than by decreasing k off .

Introduction

Since their introduction as clinical agents in the 1980s, monoclonal antibodies (mAbs) have become one of the fastest growing classes of therapeutic modalities with 24 mAbs currently on the market in the US [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] and more than 200 in clinical trials. The increasing interest in biologicals in general, and mAbs in particular, is partly due to their high success rate during preclinical and clinical development [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF]. It has been suggested that an important element that has contributed to this success is the fact that, compared to small molecules, the pharmacokinetic-pharmacodynamic (PKPD) properties of mAbs are more amenable to quantitative, mechanistic modelling and simulation (M&S)-based translation approaches across preclinical and clinical research [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF][START_REF] Agoram | Use of pharmacokinetic/pharmacodynamic modelling for starting dose selection in first-in-human trials of high-risk biologics[END_REF][START_REF] Putnam | Use of quantitative pharmacology in the development of HAE1, a high-affinity anti-IgE monoclonal antibody[END_REF][START_REF] Tabrizi | Application of quantitative pharmacology in development of therapeutic monoclonal antibodies[END_REF]. A specific example of this is the implementation of the principles of target-mediated drug disposition (TMDD) into quantitative pharmacological models to describe and predict PKPD behaviour of mAbs. Levy [START_REF] Levy | Pharmacologic target-mediated drug disposition[END_REF] introduced the term TMDD to describe the observations that for certain potent and selective drugs the fraction (relative to the dose) bound to the pharmacological target may be so high that it influences their disposition, i.e. PK and PD become interdependent. Mager and Jusko [START_REF] Mager | General pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF] were the first to propose a general PK model for drugs exhibiting TMDD, which has provided the basis for extensive further studies and development (see for example [START_REF] Tabrizi | Application of quantitative pharmacology in development of therapeutic monoclonal antibodies[END_REF][START_REF] Mager | Target-mediated drug disposition and dynamics[END_REF][START_REF] Gibiansky | Target-mediated drug disposition model: approximations, identifiability of model parameters and applications to the population pharmacokineticpharmacodynamic modeling of biologics[END_REF][START_REF] Peletier | Dynamics of target-mediated drug disposition[END_REF][START_REF] Davda | Properties of a general PK/PD model of antibody-ligand interactions for therapeutic antibodies that bind to soluble endogenous targets[END_REF]). Experimental confirmation of TMDD model predictions was provided in a recent elegant study by Abraham et al. [START_REF] Abraham | Type I interferon receptor is a primary regulator of target-mediated drug disposition of interferon-beta in mice[END_REF], who showed marked differences in interferon (INF)-β PK in wild-type and type-1 INF α/β receptor knockout mice.

Although originally proposed to describe the effects of extensive drug-target binding in tissues, TMDD has received most interest as a saturable clearance mechanism (due to, for example, receptor-mediated endocytosis; see [START_REF] Krippendorff | Nonlinear pharmacokinetics of therapeutic proteins resulting from receptor mediated endocytosis[END_REF]) for biologics, specifically mAbs. One specific example of the impact of mechanism-based PKPD models based on TMDD principles in drug discovery and development comes from the area of anti-IgE mAbs. The first generation molecule, omalizumab (Xolair R [START_REF] Hayashi | A mechanism-based binding model for the population pharmacokinetics and pharmacodynamics of omalizumab[END_REF]) is used for the treatment of allergic asthma and there has been an interest to develop more potent, second-generation mAbs with improved clinical efficacy profile [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF][START_REF] Putnam | Use of quantitative pharmacology in the development of HAE1, a high-affinity anti-IgE monoclonal antibody[END_REF]. Specifically, TMDD models have been used to explore the impact of mAb affinity for IgE in relation to in vivo potency (the dose required for a given clinical effect). Through a sensitivity analysis, Agoram et al. [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF] predicted that a ten-fold increase in mAb affinity for IgE would result in an approximately two-fold reduction in dose compared to omalizumab but that further increases in affinity were not predicted to result in additional potency improvements. In contrast, simulations did suggest that increased on-rate of mAb binding to IgE could have additional impact on the in vivo potency [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF]. However no formal analysis was presented to underwrite this hypothesis. Similarly, Sarkar et al. [START_REF] Sarkar | Cell-level pharmacokinetic model of granulocyte colonystimulating factor: implications for ligand lifetime and potency in vivo[END_REF] demonstrated the use of mechanistic models for the design of optimal biologic therapeutics such as GCSF. Therefore, some examples exist in the literature on the exploration of the affinity-potency relationships using detailed mathematical models of cellular processes, but these analyses are heuristic in nature. A systematic analysis of how mechanistic TMDD models can be employed in potency estimation of candidates is not yet available. Inspired by the beforementioned case studies and recent examples of how a rigorous mathematical analysis can aid our understanding of complex pharmacological systems and provide tools to predict essential PKPD properties of mAbs [START_REF] Mager | Quasi-equilibrium pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF][START_REF] Peletier | Dynamics of target-mediated drug disposition[END_REF][START_REF] Krippendorff | Nonlinear pharmacokinetics of therapeutic proteins resulting from receptor mediated endocytosis[END_REF] we explore the behaviour of a TMDD model with respect to the relationship between the target affinity of the mAb and its in vivo potency.

Section 1 gives the background for this study and Section 2 describes the model development from the reaction mechanism. The mathematical analysis of the resulting kinetic rate equations is depicted in this section for easy manipulation and understanding of the model. In Sections 3 and 4, the results and discussion of the model dynamics are presented and mathematical expressions to approximate the potency of the drug are derived. Section 5 gives a further analysis and validation of the relationships derived in the previous sections. Finally, the conclusions and recommendations arising from these studies are given in Section 6.

Model Equations

In this study, a one-compartment model based on the original work of Levy [START_REF] Levy | Pharmacologic target-mediated drug disposition[END_REF] will be employed in which the ligand L (drug) binds reversibly with the receptor R to form a receptor-ligand complex P as shown in the scheme in Fig. 1. The TMDD model assumes a mechanism-based reaction to explain the drug-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. 1 by 'In').

From the mechanism of the TMDD reaction shown in Fig. 1 the mathematical model can be derived using the Law of Mass Action giving the differential equations

dL dt = -k e(L) L -k on LR + k off P ( 1
)
dR dt = k in -k out R -k on LR + k off P ( 2
)
dP dt = k on LR -k off P -k e(P) P (3)
The steady state of this system is given by L = P = 0, R = k in /k out . Adding the bolus injection gives the initial conditions

L(0) = L 0 , R(0) = R 0 = k in k out , P(0) = 0. (4)
We non-dimensionalise these equations by defining the dimensionless variables

x = L L 0 , y = R R 0 , z = P R 0 , τ = k out t.
We also define the non-dimensional parameter

μ = R 0 L 0 . ( 5
)
This non-dimensionalisation is different from the one used by Peletier and Gabrielsson [START_REF] Peletier | Dynamics of target-mediated drug disposition[END_REF] who defined the new time variable as τ = k on R 0 t . The reason for this is that later on we want to explore the limits of k on → 0 and k on → ∞, and so we do not want k on to be used in the rescaling of time. We also note that k out must be non-zero for the steady state R 0 to exist, and so this is an obvious choice of parameter to use for the non-dimensionalisation. In terms of these non-dimensional quantities, equations (1)-(3) become

dx dτ = -K 1 x -K 2 xy + μK 3 z ( 6
)
dy dτ = 1 -y - K 2 μ xy + K 3 z ( 7
)
dz dτ = K 2 μ xy -(K 3 + K 4 ) z ( 8
)
with initial conditions

x(0) = 1, y(0) = 1, z(0) = 0,
where the dimensionless parameters are defined by

K 1 = k e(L) k out , K 2 = k on R 0 k out , K 3 = k off k out , K 4 = k e(P) k out .
The constant term in [START_REF] Mager | General pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF] is derived by using the definition of the steady state R 0 in (4). We note that the six parameters and one initial value (L 0 ) of the dimensional equations have been reduced [START_REF] Sun | Mechanism-based pharmacokinetic/pharmacodynamic (PK/PD) modelling for biotechnology products, poster presented at Advanced Methods of PKPD Systems Analysis[END_REF].

to five non-dimensional constants for the non-dimensional equations, which simplifies the model equations while keeping the essence of the system.

We use parameter values from the IgE mAb omalizumab case study [START_REF] Sun | Mechanism-based pharmacokinetic/pharmacodynamic (PK/PD) modelling for biotechnology products, poster presented at Advanced Methods of PKPD Systems Analysis[END_REF] (which were also used in [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF]) as an example. The numerical values of the dimensioned and dimensionless parameters for this case are given in Table 1, where we have also included the value of the equilibrium dissociation constant K D = k off /k on . For the dimensioned parameters, we use these units hereafter without mentioning them specifically. Using these parameters, the time profile of the ligand, receptor and ligand-receptor complex are shown in Fig. 2 for a short time (just past the minimum of R ) and for a longer time interval in Fig. 3. Clearly, plots of the non-dimensional quantities x, y and z will be similar but with a difference scale on the vertical axis.

These numerical results show different phases occurring. Initially, in a very short time, the receptor R drops down to a low value while the product P shows a corresponding jump. There is also an initial sharp but small drop in the ligand L. This phase is dominated by the binding action of the ligand to the receptor resulting in the product. In the next phase, which happens on a much slower timescale, the receptor and the product gradually increase, while the ligand continues with a gradual decrease. Finally, the system settles back to its steady state values.

Approximation of the Drug Potency -Method 1

In our analysis, the in vivo potency of the drug is defined to be the minimum free receptor level that can be achieved for a particular dose of the drug. Potency is defined as the amount of drug (the dose) that is needed to produce a defined effect [START_REF] Gabrielsson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF]. We use the minimum free receptor level as the measure of 'defined effect' since this can conceptually be related to a level of clinical outcome. For example, in the case of omaluzimab and follow-on IgE mAbs, reducing mean free IgE levels to approximately 10 IU/ml in asthma patients has been used as a target to achieve the desired clinical efficacy [START_REF] Putnam | Use of quantitative pharmacology in the development of HAE1, a high-affinity anti-IgE monoclonal antibody[END_REF].

Mathematically, this minimum receptor level, R min , is obtained by finding the minimum point on the receptor-time profile (see Fig. 2(b)) and is therefore defined as the point where the time derivative of R (or y in dimensionless form) vanishes. Setting ẏ = 0 in [START_REF] Mager | General pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF] implies that this minimum occurs when 1y -

K 2 μ xy + K 3 z = 0, which gives y = μ(1 + K 3 z) μ + K 2 x .
The problem with this is that we do not know the values of x and z when the minimum occurs, and so all this gives us is a two-dimensional surface in the three-dimensional phase space on which the minimum must occur but the position on this surface where the actual minimum occurs for the given trajectory is not known. To address this problem, we rewrite equations ( 6)-( 8) in terms of different variables. In particular, we introduce a non-dimensional form of the total amount of ligand (free and bound) together with the total amount of receptor (free and bound). Thus, we define

u = L + P L 0 = x + μz v = R + P R 0 = y + z
We now rewrite our equations in terms of the variables u , v and y , giving

u = -K 1 u -μ(K 4 -K 1 )(v -y) ( 9 ) v = 1 -y -K 4 (v -y) (10) ẏ = 1 -y - K 2 μ uy + (K 2 y + K 3 )(v -y) ( 11
)
with initial conditions

u(0) = 1, v(0) = 1, y(0) = 1.
The concentration-time profile for the total ligand L + P and the total receptor R + P are shown in Fig. 4. We assume that the factor μ/K 2 = k out /(k on L 0 ) is small which requires that k on L 0 is large relative to k out . This condition will be satisfied if the amount L 0 of ligand injected is sufficiently large or if the ligand and the receptor combine very fast, i.e., k on is large. We define

ε = μ/K 2 ,
as our small parameter.

As before, we can set ẏ = 0 in (11) which, using our new parameter ε , gives

1 -y - 1 ε uy + (K 2 y + K 3 )(v -y) = 0.
We note that in this case, we have a quadratic equation for y which has one positive and one negative solution. The positive solution is given by

y = - 1 2εK 2 [u + ε(1 + K 3 -K 2 v)] -[(u + ε(1 + K 3 -K 2 v)) 2 + 4ε 2 K 2 (1 + K 3 v] .
Since we have assumed that ε is small, the first term under the square root will dominate the second. Thus, taking this first term out of the square root and then expanding the square root term gives

y = ε 1 + K 3 v u + ε(1 + K 3 -K 2 v) + O ε 2 . ( 12
)
Of course, we still have the problem that this formula for y involves the variables u and v . However, the advantage of working with equations ( 9)-( 11) is that the variables u and v evolve on a much slower timescale than y . To see this more clearly, we introduce the fast time variable τ = t/ε . With this time variable, the system (9)-( 11) becomes

u = -ε [K 1 u + εK 2 (K 4 -K 1 )(v -y)] (13) v = ε [1 -y -K 4 (v -y)] ( 14
)
y = -uy + ε [1 -y + (K 2 y + K 3 )(v -y)] , ( 15
)
where u stands for du/dτ , etc. We note that

u (0) = -εK 1 , v (0) = 0, y (0) = -1.
Provided that εK 1 = k e(L) /(k on L 0 ) is small, which is the case for the parameter values we are considering, then the derivative of y at time zero is much greater than the corresponding derivatives of u and v , and hence y will initially change at a much faster rate than u and v , as can be seen from Figs 2 and4. Indeed, the variables in Fig. 4 do not show the initial sharp change that can be observed for R (or equivalently y ) in Fig. 2. Provided that εK 4 = k e(P) /(k on L 0 ) and ε = k out /(k on L 0 ) are also small, the derivatives u and v are small as they have an overall factor ε . So during an initial fast phase, u and v will remain approximately constant while y quickly decreases to its minimum and so we can approximate this minimum by setting u = u(0) = 1 and v = v(0) = 1. Substituting these values into (12) and ignoring the higher order terms gives an approximation of the minimum value of y as

y (1) min = ε 1 + K 3 1 + ε(1 + K 3 -K 2 ) . ( 16
)
Converting this expression back to the dimensional variables, we obtain an approximation for the minimum of R as

R (1) min = R 0 (k off + k out ) k on (L 0 -R 0 ) + k off + k out . ( 17
)
The justification of this approximation is obtained as a consequence of the reaction between L and R , which produces the product P , happening on a much faster timescale than the other dynamic processes. From the non-dimensional equations ( 13)-( 15), we can see that the reduction of y will happen on a fast timescale if ε is small and u and v are not depleted too quickly by elimination relative to the reaction rate. This requires that k e(L) , k e(P) and k out are small relative to k on L 0 . Clearly, these conditions can be fulfilled in two ways, namely by L 0 being large relative to k out , k e(P) and k e(L) or by k on being large relative to k out , k e(P) and k e(L) (with no corresponding requirement that L 0 also be large).

Looking at the limiting cases, we see that varying k on gives

lim kon→0 R (1) min = R 0 , lim kon→∞ R (1) min = 0. ( 18
)
Similarly, varying k off gives lim

koff →∞ R (1) min = R 0 , lim koff →0 R (1) min = R 0 k out k on (L 0 -R 0 ) + k out , ( 19
)
We note that in the limit as K D = k off /k on → ∞ there is no reaction between the ligand and the receptor and hence no product is formed. Thus, R remains at its steady state value R 0 and so we have obtained the correct limit in this case.

In the alternative case when K D = k off /k on → 0, we note that increasing k on results in the minimum level of the receptor decreasing to zero, while decreasing k off results in the minimum receptor level decreasing to a non-zero level.

We note in passing that the formula for R

(1

)
min can be expressed in terms of the single parameter (k off + k out )/k on for given L 0 and R 0 and so an increase in k on will have the same effect as a corresponding decrease in k off + k out . However, k out is not a constant over which we have any control, and so this observation is not particularly helpful in practice.

Approximation of the Drug Potency -Method 2

In the previous section, we worked with the dimensionless variables for the total ligand (u), total receptor (v) and free receptor (y) . We now consider a similar approach using instead the dimensionless variables for the total ligand (u) , free receptor (y) and product (z). In this case, the condition for the minimum of y (i.e. ẏ = 0) gives

ε(1 -y) -(u -εK 2 z)y + εK 3 z = 0,
which is linear in y . Solving this equation for y we obtain

y = ε(1 + K 3 z) u + ε(1 -K 2 z) .
As before, we note that u remains approximately constant during the initial fast phase and so we again take u = 1 . However, z changes rapidly during the initial phase, as does y , and so we cannot use the initial value of z . The assumption used previously that v = 1 implies that y+z = 1 or z = 1y . Now at the minimum point, if we assume that y = O(ε), then z = 1 -O(ε) at this point. Since ε is assumed to be small, we simply take z = 1 as the leading order approximation. With these two assumptions, we obtain an approximation for the minimal receptor level as

y (2) min = ε 1 + K 3 1 + ε(1 -K 2 ) , ( 20
)
which is very similar, but not quite the same, as the previous approximation y

(1)
min given in [START_REF] Hayashi | A mechanism-based binding model for the population pharmacokinetics and pharmacodynamics of omalizumab[END_REF]. Converting back to the original coordinates gives

R (2) min = R 0 (k out + k off ) k on (L 0 -R 0 ) + k out . ( 21
)
In the derivation of this approximation, we used the same assumptions as for R

(1) min , i.e., k e(L) , k e(P) and k out are small compared to k on L 0 , plus the extra assumption that k off is small compared to k on L 0 .

It is easily verified that y

(

min + O(ε 2 ). We also note that

lim kon→∞ R (2) min = lim kon→∞ R (1) min = 0, (22)
and that lim

koff →0 R (2) min = lim koff →0 R (1) min = R 0 k out k on (L 0 -R 0 ) + k out , ( 23
)
and so the limiting values for the two approximations R

(1) min and R

(2) min as K D = k off /k on → 0 are the same. However, in the other limit with

K D = k off /k on → ∞, we find that lim kon→0 R (2) min = R 0 1 + k off k out , lim koff →∞ R (2) min = ∞.
Thus, the correct values are not obtained in this limit. This is expected since we assumed above that y min = O(ε), which in turn implies that k off /(k on L 0 ) is small. When k off /k on is large, this condition is violated and the approximation is not valid. We note that even for a large (fixed) value of L 0 , if K D becomes sufficiently large, then k on L 0 will no longer be large relative to k off and so this approximation fails.

Validation and Further Analysis of the Approximations

We now compare these predicted approximate values with numerical values obtained from the differential equations. We use the values of the constants from the omalizumab case study, as given in Table 1. Since ε = 0.0938 is quite small and K D = 1.52, both approximations will be valid in this case. Solving the differential equations numerically gives the true minimum of R to be R min = 0.5811. Our two approximations give values of K D = 1.52. To make it easier to compare results for varying either k on or k off , we first fix k on = 0.592 and vary k off via k off = 0.592 K D and then we fix k off = 0.9 and vary k on = 0.9/K D . The resulting R min data points from the numerical simulation of the differential equations and approximation curves R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min and R [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF] min are plotted as functions of the variable K D in Fig. 5. This illustrates that R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min is a good approximation for all values of K D . For small values of K D , the approximation R [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF] min is better than R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min , but for larger values of K D , the approximation R (2) min = 0.2765, while the numerical simulation gives R min = 0.3014. For k off = 0.9 and k on → ∞, we have both the simulation and the approximations converging to zero.

It can be seen that the approximation R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min corresponds well to the results obtained by simulation over the whole range of K D values, while the approximation R [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF] min is a good approximation when R min /R 0 is small, which occurs when either L 0 or k on is large. Thus the analytical expressions capture the relationship between the potency ( R min ) and the affinity constants (k off and k on ) very well.

Another way to view the expression for R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min is by plotting it as a function of the two variables log 10 (k off ) and log 10 (k on ) as shown in Fig. 6. The R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min contours of this surface are shown in Fig. 7, while the contours with k off and k on constant are shown in Figs 8 and9 respectively. We note in Fig. 7 that when log 10 (k off ) is positive, the contours of R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min appear to be approximately straight lines of slope one. On the other hand, when log 10 (k off ) is less than -1, the contours are approximately horizontal lines. This can be explained as follows.

If k off k out , then k off will make a negligible contribution to the numerator and the denominator of R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min and so we can set it to zero, giving

R (1) min ≈ R 0 k off k on (L 0 -R 0 ) + k off = R 0 K D L 0 -R 0 + K D , for k off k out .
Thus, in this case, R Taking logs, we then find that

K D = R c (L 0 -R 0 ) R 0 -R c , for k off k out .
log 10 (k on ) = log 10 (k off ) -log R c (L 0 -R 0 ) R 0 -R c , for k off k out
and so the relationship between log 10 (k off ) and log 10 (k on ) is linear with slope one, as seen in Fig. 7, for k off sufficiently large. At the other extreme, when k off k out , then k off will be negligible compared to k out and so, setting it to zero, we obtain

R (1) min ≈ R 0 k out k on (L 0 -R 0 ) + k out , for k off k out .
In this case, setting R

(1) min = R c and solving for k on , we obtain

k on = k out (R 0 -R c ) R c (L 0 -R 0 ) , for k off k out
and clearly this expression for k on does not depend on k off (since we set it to zero!) and so log 10 (k on ) is a constant. Thus, the contours asymptote to a constant as log 10 (k off ) → -∞ (k off → 0) and asymptote to a straight line of slope one as log 10 (k off ) → ∞ ( k off → ∞). The in-between region, when k off and k out are of similar order, gives the curve that joins these two straight lines.

We are now able to derive further information on the effect of varying k off or k on on R min . When k off is large relative to k out , we have seen that R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min effectively depends only on the single parameter K D , in which case increasing k on by a factor of α will have the same effect on R min as decreasing k off by a factor of α . However, when k off is either small or of similar magnitude to k out (as is the case for omalizumab), then increasing k on by a factor of α will have a much greater effect than reducing k off by a factor of α .

We note that over the whole range of values of log 10 (k off ) to move from one contour in Fig. 7 to a lower one can always be achieved by a relatively small increase in log 10 (k on ) . However, moving to a lower contour by decreasing log 10 (k off ) can easily be achieved with a relatively small decrease when log 10 (k off ) is large, but requires a much larger decrease for smaller values of log 10 (k off ), and for smaller values still, it is not possible to move to a lower contour by decreasing log 10 (k off ). This is consistent with the results shown in Fig. 5.

Finally, we consider what increase in k on , or decrease in k off , is required to reduce R min by a factor of two (which correponds to increasing the efficacy from say 90% to 95%). We suppose that R

R c = R 0 (k off + k out ) k 0 on (L 0 -R 0 ) + k off + k out , R c 2 = R 0 (k off + k out ) k 1 on (L 0 -R 0 ) + k off + k out .
Eliminating R c from these equations and solving for k 1 on gives

k 1 on = 2k 0 on + k off + k out L 0 -R 0 , ( 24
)
and so we conclude that to reduce R min by a factor of two, k on has to be more than doubled, since the second term on the right of (24) is positive.

If we define K 0 D = k off /k 0 on and K 1 D = k off /k 1 on , then from (24) we find that

K 1 D = 1 2 K 0 D 2k 0 on (L 0 -R 0 ) 2k 0 on (L 0 -R 0 ) + k off + k out .
Since the term in the brackets is less than one, clearly K D must be reduced by a factor greater than two in order to reduce R min by a factor or two, which is consistent with the above statement regarding k on .

A similar calculation where k off is varied rather than k on gives

k 1 off = 1 2 k 0 off - 1 2 k out - (k 0 off + k out ) 2 4k on (L 0 -R 0 ) + 2(k 0 off + k out ) . ( 25
)
In this case, to reduce R min by a factor of two, k off must be reduced by more than a factor of two. Moreover, it is quite possible that this formula gives k 1 off < 0, which of course implies in such a case that it is not possible to reduce R min by a factor of two by reducing k off , which we have already noted above.

If we define K0 D = k 0 off /k on and K1 D = k 1 off /k on , then from (25) we have that

K1 D = 1 2 K0 D - k out 2k on - (k 0 off + k out ) 2 k on [4k on (L 0 -R 0 ) + 2(k 0 off + k out )] = 1 2 K0 D - k 2 - ( K0 D + k) 2 4(L 0 -R 0 ) + 2( K0 D + k)
where k = k out /k on . The comments above regarding k off apply also to KD .

Discussion and Conclusion

The main objective of the work presented in this paper was to explore the TMDD model through a rigorous mathematical analysis with regards to the relationship between the target affinity of a mAb versus its in vivo potency. In our experience, this topic invariably gets raised in mAb drug discovery and development programs, mainly because the maximum dose for routine clinical use in patients is typically more stringently constrained for mAbs compared to small molecules due to non-oral route of administration, formulation complexities and cost of goods, and because increasing the affinity of a mAb is a time-consuming process.

The first, perhaps obvious, conclusion of our analysis is that for a given dose (L 0 ) the minimum value of free target (R min ) can be decreased (equivalent to increasing in vivo potency) by increasing k on or by decreasing k off . While this conclusion may be obvious, what may not be so apparent is that there is a significant difference in the minimum receptor obtained when K D is reduced either by reducing k off or by increasing k on . We note from [START_REF] Mager | Quasi-equilibrium pharmacokinetic model for drugs exhibiting target-mediated drug disposition[END_REF] that as k on → ∞, then R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min → 0, whereas from [START_REF] Dimitrov | Engineered CH 2 domains (nanoantibodies)[END_REF], as k off → 0, then R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min tends to a non-zero value, with the same limits being obtained from R [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF] min , as already noted. Thus, we have a saturation effect when decreasing k off , in that further reductions in the value of k off will yield only limited reductions in R min , with a positive limiting value as k off → 0. There is no such saturation when increasing k on , so that R min is consistently reduced towards the limiting value of zero for increasing values of k on .

We see from [START_REF] Dimitrov | Engineered CH 2 domains (nanoantibodies)[END_REF] and ( 23) that the limiting value of R min as k off → 0 in both cases is given by

lim koff →0 R (1) min = lim koff →0 R (2) min = k in k out k on (k out L 0 -k in ) + k 2 out ,
where we have substituted R 0 = k in /k out . Thus, we can see that this saturation level will be reduced by increasing k on or L 0 . Theoretically, the saturation level could also be reduced by decreasing k in or by increasing k out .

This sheds new light on our previous work [START_REF] Agoram | The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics[END_REF] in which, through simulations, we predicted that the maximum increase in potency that could be achieved with an anti-IgE antibody by only reducing k off would be approximately two-fold compared to omalizumab, and that this would be achieved with a five-to-ten-fold increase in affinity. Clinical data reported subsequently by Putnam et al. [START_REF] Putnam | Use of quantitative pharmacology in the development of HAE1, a high-affinity anti-IgE monoclonal antibody[END_REF] on the high-affinity anti-IgE mAb, HAE1, were consistent with this prediction, since it was shown that HAE1 achieved an approximately two-fold improvement in in vivo potency compared to omalizumab, whereas it displayed a more than twenty-fold higher affinity for IgE. Interestingly, the data presented in [START_REF] Putnam | Use of quantitative pharmacology in the development of HAE1, a high-affinity anti-IgE monoclonal antibody[END_REF] show that the affinity improvement of HAE1 compared to omalizumab was entirely driven by a reduction in k off , consistent with our prediction that this parameter is associated with a saturation effect regarding its impact on in vivo potency.

As far as we are aware there is no experimental data to support our hypothesis that such a saturation effect does not exist for k on , but our analysis does suggest that a mAb potency optimisation strategy focused on increasing k on rather than decreasing k off could be advantageous. Currently, the majority of marketed mAbs and those in clinical development are IgE's of about 150 kDa size and their k on values are generally uniform and limited by their size. However, even if rational optimisation of k on is currently not experimentally possible, our analysis indicates that between two antibodies of sufficiently low k off (in the saturation effect range), one antibody with a higher k on is more potent even if it may have higher k off and hence, potentially higher overall K D . While this conclusion may already have been derived in an empirical manner within certain areas of biologics discovery (especially with the extensive research into highly labile targets such as interleukins) we have provided the first systematic quantitative analysis that can guide rational optimisation of mAbs within the context of the TMDD framework. For example, in recent years, significant efforts have been put into the development of novel human and non-human scaffolds ('nanobodies') of much smaller size (see [START_REF] Dimitrov | Engineered CH 2 domains (nanoantibodies)[END_REF][START_REF] Van Bockstaele | The development of nanobodies for therapeutic applications[END_REF]). Although currently these efforts appear to be mainly motivated by predicted improvement of tissue penetration, systemic stability and preferential cleft recognition [START_REF] De Genst | Molecular basis for the preferential cleft recognition by dromedary heavy-chain antibodies[END_REF], an intriguing question that follows on from the present analysis is whether nanobodies could display faster k on rates due to their smaller size and therefore may be more amenable to optimisation of in vivo potency compared to traditional mAbs.

We have chosen our measure of drug potency to be the minimum free receptor level, although there are other measures that could be used, such as the maximum of the complex P . However, this maximum typically occurs at a much later time than the minimum of R (at approximately 4.5 days for the example shown in Fig. 3(b)) and at this later time, our assumptions that u = v = 1 are no longer valid. Thus, the methods that we have described are not applicable in this case.

Alternative measures of potency that could be considered are (i) the maximum of the receptor occupancy P/R or (ii) the minimum of the proportion of free receptor relative to the total amount of receptor R/(R + P ). In terms of the non-dimensional variables, these quantities are z/y and y/v respectively. To find the maximum/minimum levels, we again differentiate and set to zero. For the receptor occupancy, we then obtain Thus, the equation to be solved is the same in both cases. Reworking method 1 above for this equation, we find that the maximum/minimum occurs when

R = R 0 (k out + k off ) k on (L 0 -R 0 ) + 2k out + k off -k e(P)
Clearly, this is very similar to our formula for R [START_REF] Nelson | Development trends for human monoclonal antibody therapeutics[END_REF] min given by [START_REF] Sarkar | Cell-level pharmacokinetic model of granulocyte colonystimulating factor: implications for ligand lifetime and potency in vivo[END_REF]. In particular, the broad conclusions that there is saturation to a non-zero value as k off → 0 and no such saturation as k on → ∞ still hold in both cases.

It should be pointed out that an important assumption in our analysis is that k on L 0 must be significantly larger than k e(L) , k e(P) and k out . It is easily verified that the quotients k e(L) /(k on L 0 ), k e(P) /(k on L 0 ) and k out /(k on L 0 ) are all less than 0.1 for the IgE mAb omalizumab case study and hence are sufficiently small for the analysis to work well. However, while this assumption may generally hold true for blocking and neutralising antibodies, it may not be true for agonists and in cases where pharmacological effect may be exerted at low levels of target occupancy. Therefore, our conclusions may not be valid in those cases. However, we do not believe this greatly limits applicability of the simplification we have provided, since in the cases where the underlying assumption is violated, lower doses are only required for clinical efficacy and therefore affinity requirements tend to be less stringent. It should also be noted that this analysis is only applicable for a rather simple antigen-antibody system which interacts without diffusion barriers, avidity effects, and other complexities. These complicating factors should be considered to be able to extend the model to more realistic drug discovery situations.

In conclusion, our analysis of the TMDD model has provided a mathematical framework that relates intrinsic pharmacological and pharmacokinetic properties of mAbs to their in vivo potency. The finding that the greatest potency improvements can be achieved, at least in theory, through modulation of k on could provide a basis for new strategies to drive the discovery of a new generation of mAbs. The simple formula we have provided can also substantially reduce the need for complex PKPD analyses resources to get an initial estimate of required affinity at the lead optimisation stage.

Mathematical Analysis of the Pharmacokinetic/Pharmacodynamic (PKPD)

Behaviour of Monoclonal Antibodies: Predicting in vivo Potency P.J. Aston, G. Derks, A. Raji, B.M. Agoram, P.H. van der Graaf

Highlights

x Potency is measured by the minimum free receptor following an injection of ligand.

x Two approximations of this minimum are derived from the differential equations.

x A saturation effect is demonstrated when the dissociation constant k off is reduced.

x There is no such saturation when the association constant k on is increased.

Figure 1 :

 1 Figure 1: The TMDD reaction mechanism

Figure 2 :

 2 Figure 2: Concentration-time profile for 0 ≤ t ≤ 1 of (a) the free ligand L , and (b) the free receptor R (blue) and the ligand-receptor complex P (red).

Figure 3 :

 3 Figure 3: Concentration-time profile for 0 ≤ t ≤ 20 of (a) the free ligand L , and (b) the free receptor R (blue) and the ligand-receptor complex P (red).

Figure 4 :

 4 Figure 4: Concentration-time profile of (a) the total ligand L + P and (b) the total receptor R + P . Note the small variation in values on the vertical L + P -axis.

Figure 5 :

 5 Figure 5: The minimal value of R as a function of K D with (a) K D plotted on a log scale for a wide range of K D values and (b) K D plotted on a linear scale for K D small. In both plots, the red solid curve ( R (1) min ) and dash-dotted curve ( R (2) min ) and points * (numerical simulation) are for k off = 0.9 and kon = 0.9/K D . Similarly, the blue solid curve ( R (1) min ) and dashed curve ( R (2) min ) and points • (numerical simulation) are for kon = 0.592 and k off = 0.592 K D .

( 2 )

 2 min starts to fail as expected. For k off = 0 and k on = 0.592, we have the approximation R

( 1 )

 1 min depends only on the single parameter K D = k off /k on and not separately on k off and k on . If we now consider a contour which corresponds to a constant value of R min = R c , for some given constant R c , then we can solve the equation R (1) min = R c for K D giving

Figure 6 :

 6 Figure 6: Plot of R (1) min given by (17) as a function of log 10 (k off ) and log 10 (kon) .

Figure 7 :

 7 Figure 7: Contours of the surface shown in Fig. 6 for constant values of R min .

Figure 8 :

 8 Figure 8: Contours of the surface shown in Fig. 6 for constant values of k off .

Figure 9 :

 9 Figure 9: Contours of the surface shown in Fig. 6 for constant values of kon .

( 1 )

 1 min = R c when k on = k 0 on and that R (1) min = R c /2 when k on = k 1 on . This gives the two equations

Table 1 :

 1 Numerical values of the dimensioned and dimensionless constants for IgE mAb omalizumab

	Kinetic Constant	Value	Parameter	Value
	k e(L) k e(P) k out k off k on	0.024 day -1 0.201 day -1 0.823 day -1 0.900 day -1 0.592 (nM day) -1	K 1 K 2 K 3 K 4 μ	0.02916 1.93353 1.09356 0.24423 0.18144
	K D k in	1.520 nM k out R 0 nM day -1
	L 0	14.8148 nM
	R 0	2.688 nM

Acknowledgements

This work was financially supported through the Biopharma Skills Project of the Universities of Surrey and Reading, jointly funded by the Higher Education Funding Council for England's Economic Challenge Investment Fund (ECIF) and the South East England Development Agency (SEEDA).

We are grateful to Prof Bert Peletier (Leiden) for helpful discussions and encouragement. We also thank Dr David Lloyd (Surrey) for his interest and advice. We thank the referee for some helpful comments that have improved the paper.