Do we really understand how drug eluted from stents modulates arterial healing?

of drug on SMC growth in vitro . Our results highlight that, at least for Sodium Salicylate and Paclitaxel, the current state-of-the-art nonlinear saturable binding model is incapable of capturing the proliferative response of SMCs across a range of drug doses and exposure times. Our findings potentially have important implications on the interpretation of current computational models and their future use to optimise and control drug release from DES and drug-coated balloons.


Evolution of drug-eluting stents
Drug-eluting stents (DES) were introduced over two decades ago to counteract the high in-stent restenosis (ISR) rates associated with their bare-metal (BMS) predecessors (Bell et al., 1992;Landau et al., 1994;van der Hoeven et al., 2005). Initial designs incorporated either sirolimus (SIR) or paclitaxel (PTX) within a durable polymer that was coated onto a stainless steel platform. The clinical results with these initial DES were remarkable, with the rate of ISR reducing from the 17-41% associated with bare metal stents to below 10% (Moses et al., 2003;Park et al., 2003;Buccheri et al., 2016). Amid concerns regarding the increased stent thrombosis (ST) risk (2.2% for DES and 1.5% for BMS (Tada et al., 2013) and delayed arterial healing associated with these initial DES designs (Byrne et al., 2014;Nakazawa et al., 2008;Byrne et al., 2009;Holmes et al., 2010), further advancements in technology focused on improved biocompatibility and thinner struts (Akin et al., 2011;Puranik et al., 2013). This second generation of DES resulted in a modest improvement in thrombosis and restenosis rates of 1% and 6%, respectively (Byrne et al., 2014;Whitbeck and Applegate, 2013).
Further generations have featured alternative drugs, predominantly analogues of sirolimus (e.g. everolimus, zotarolimus, biolimus) and biodegradable polymer coatings, generally resulting in more favourable results when compared to PTX-eluting stents, yet non-inferior when compared to sirolimus (Stone et al., 2010;Navarese et al., 2014;Byrne et al., 2011). Most recently, polymer-free stents and fully resorbable stents have been developed, with the Biofreedom stent an example of the former, a superior alternative to BMS in the context of high-risk bleeding patients requiring dual anti-platelet therapy (DAPT) (Urban et al., 2015), whilst the Absorb stent (an example of the latter) has been removed from the market, amid ongoing safety concerns related to increased rates of thrombosis (Cassese et al., 2016).
It is indisputable that DES work well in the majority of cases. However, despite all the aforementioned technological innovations, complications continue to occur, with repeat revascularisation procedures still required in approximately 6% of cases within 2 years (Byrne et al., 2015;Byrne et al., 2018), and even more so in cases associated with complex lesions (Serruys et al., 2009).
Delayed healing has been a recurring problem in DES, where durable coatings were associated with hypersensitivity, impacting the reendothelization process (Joner et al., 2006). Various efforts have presented alternate results, where the occurence of adverse events and neointima development is comparable between durable and resorbable coated DES in some studies (Byrne et al., 2011;Yuan et al., 2014;von Birgelen et al., 2016;Zocca et al., 2018;Guagliumi et al., 2018;Kereiakes et al., 2019) yet inferior in others with respect to late-stent thrombosis and early strut coverage (Kobayashi et al., 2017;Nojima et al., 2019;Serruys et al., 2013;Stefanini et al., 2011). Conflicting evidence considering various DES, supported by a recent meta-analysis (Taglieri et al., 2020), suggests that the drug and release kinetics may be instrumental in cases of delayed arterial healing and late stent thrombosis (Nakazawa et al., 2008).

Mathematical and computational modelling in DES research
Over the past two decades, mathematical and computational modelling has emerged as a powerful tool to help understand the performance of DES. Such research has predominantly focused on three key areas: (i) structural mechanics and computational fluid dynamics (CFD); (ii) drug elution kinetics, uptake and retention in tissue, and (iii) the arterial healing response.
The most advanced models of type (i) are now capable of describing stent deployment and the subsequent alteration to arterial fluid dynamics in 3D patient-specific arteries. Moreover, these models have enabled correlations between haemodynamic indices, such as wall shear stress (WSS) and Oscillatory Shear Index (OSI) and clinical outcome (i.e. restenosis) (Morlacchi et al., 2011;Chiastra et al., 2013;Chiastra et al., 2016;Colombo et al., 2020).
State-of-the art models of type (ii) have elucidated the importance of drug release kinetics and retention on patient outcome (McGinty, 2014;Tzafriri et al., 2012;Rossi et al., 2012;McGinty and Pontrelli, 2016;Artzi et al., 2015;Zhu and Braatz, 2014;McKittrick et al., 2019). Whilst drug release and subsequent transport has been modelled extensively, in increasingly sophisticated modelling frameworks, the consequent effect of spatiotemporal drug levels on cell behaviour has largely been neglected, which is somewhat surprising given the key role that cells play in neointimal hyperplasia formation and the overall arterial healing response.
Lastly, a series of models of type (iii) have emerged that aim to simulate the healing response of the artery following stent deployment. The majority of these consider agent based models (ABMs) (Tahir et al., 2011;Tahir et al., 2013;Tahir et al., 2014;Tahir et al., 2015;Zun et al., 2017;Pavel Zun et al., 2019;Boyle et al., 2010), whereby a set of rules are prescribed to dictate cell behaviour. In some cases, ABMs are coupled with discrete models of flow (Tahir et al., 2011;Tahir et al., 2013;Tahir et al., 2014;Tahir et al., 2015;Zun et al., 2017;Zun et al., 2019), enabling the calculation of flow indices such as wall shear stress (WSS) to be included within the rules. The most common alternative to ABMs is continuum models typically consisting of several coupled reaction-diffusion equations that depict how a number of species (e.g. growth factors, cells, etc.) contribute to the growth of restenotic tissue, either in the presence (Boland et al., 2019) or absence of flow (Lally and Prendergast, 2006;Boland et al., 2016;Boland et al., 2017;Escuer et al., 2019;Fereidoonnezhad et al., 2017;He et al., 2020). Whilst these ABM and continuum models have been useful in gaining insight into the physical mechanisms of tissue remodelling following stent deployment, they lack the action of drug on species evolution; the key driver in preventing the re-occlusion of the vessel.

Current models of smooth muscle cell proliferation inhibited by drug
Currently, mathematical models which detail the effect of drug on cell proliferation predominantly emanate from cancer applications (Enderling and Chaplain, 2014), where the purpose of the drug and delivery method is somewhat different. In such scenarios, the drug is administered via a range of methods (orally, intravenously) routinely to kill tumour cells. In stents, the drug is applied through a sustained local delivery and is intended to induce a cytostatic response.
Often emulating in vitro environments, the existing models aim to capture how the drug halts cell proliferation. The simplest approaches utilise coupled linear ordinary differential equations (ODE's), stemming from the Skipper-Schabel models, investigating two cell types; resting and proliferating (Stein et al., 2018;Panetta, 1997;Skipper et al., 1970). Others consider a more complex approach to exhibit various effects, including: (i) the evolution of damage towards death of cells (Del Bene et al., 2009); (ii) cell-cycle progression (Sciumè et al., 2013;Benzekry et al., 2012;Roe-Dale et al., 2011); and (iii) drug resistance (Roe-Dale et al., 2011). More recent work suggests different models may be necessary to portray the effect different drugs have on cell growth (Jarrett et al., 2019). The reader is referred to recent review articles for further information on mathematical approaches aimed to applications in cancer (Enderling and Chaplain, 2014;Ribba et al., 2014;Jarrett et al., 2018).
Over two decades ago, Schwartz et al. (1996) presented a simple model to describe exponential smooth muscle cells (SMCs) growth, considering the total number of cells and the number of proliferating cells. They achieved good agreement with ex vivo data from rat carotid arteries and noted that this behaviour is highly species-dependent. However, to the best of our knowledge, very few models exist that depict the action of drug in the context of SMC growth, most of which have been based upon the large literature on cell proliferation in the context of cancer modelling, Rossi et al. (2012), as part of a more complicated model of bioresorbable stent kinetics, adopted a corpuscular approach and described the time evolution of number of cells through a population balance equation including terms to account for cell growth, replication by mitosis and death due to apoptosis and necrosis. They solved these equations using the method of moments, specifying a unimodal lognormal distribution function for cell masses and allowing cells to divide only when they reach some critical mass. Although the functional form of the growth, mitosis and death rates is not provided in the article, through a fitting procedure, Rossi et al. (2012) demonstrate reasonable agreement between their model and existing data of SMC proliferation subject to drug exposure in vitro. In line with their previous modelling of avascular tumour growth, the authors modelled the effect of drug on cell growth assuming a simple linear dependence, i.e. drug reduces the rate of growth of cells in proportion to concentration of drug via some rate constant. Rossi et al. (2012) proceeded to apply their model to the in vivo situation, simulating neointimal thickness as a function of time by relating the number of cells to the thickness of tissue they would produce through a simple equation involving the average radius of a single cell layer and the maximum number of cells allowed in each proliferating layer.
More recently, Peddle et al. (2018) attempted to adapt models of cancer biology to model the biology of restenosis and the effect of drug elution. They presented a partial differential equation (PDE) model, detailing proliferating and quiescent SMC fractions as separate species.
Growth and death within each phase and transition between phases was accounted for by linearly combining a number of rate constants. Vessel thickening was then simulated by imposing a Stefan condition. The action of the drug was described in one of two ways: either through a growth inhibition model or through a transition-blocking model.
The state-of-the-art models of drug transport within arterial tissue include specific binding of drug to target cell receptors, following the pioneering works of Tzafriri et al. (2012), where it was implicated that specific receptor saturation was a key determinant of success (i.e. reducing excessive neointimal hyperplasia), at least for SIR. Similar nonlinear saturable binding models have been proposed for PTX (Tzafriri et al., 2009). Whilst these non-linear binding models are now used routinely in computational studies investigating drug-elution and retention in tissue, there is a disconnect in that the aforementioned models of Rossi et al. (2012) and Peddle et al. (2018) assume the effect of drug on SMCs follows somewhat simpler kinetics. In this article we seek to test the ability of mathematical models of varying complexity to capture existing published experimental data of cell proliferation subject to drug delivery.

Outline
In Section 2 we discuss the role of anti-proliferative drugs on cell growth and provide an overview of existing experimental (in vitro) studies in the literature, before discussing in detail the specific data sets that we consider in this work. We then, in Section 3, present our mathematical models of cell growth in the absence and presence of drug. We discuss briefly how we solve the governing equations and outline our parameter estimation methodology. In Section 4 we provide results comparing each mathematical model to the experimental data. Finally, in Section 5 we discuss the implications of our work for both the mathematical modelling and experimental communities in the field.

Brief Review of in vitro experimental data
It has been well documented that drugs such as PTX and SIR inhibit cell growth in a dose-dependent manner in vitro (Oberhoff et al., 2002;Axel et al., 1997;Wessely et al., 2006;Blagosklonny et al., 2004;Wiskirchen et al., 2004). As illustrated by Fig. 1, these are cell-cycle specific drugs, predominantly exerting their effect at a single point in the cell cycle (Oberhoff et al., 2002;Shah and Schwartz, 2001).
PTX and SIR are both thought to enter the cell membrane through rapid diffusion (Oberhoff et al., 2002;Ndungu et al., 2010). The former binds specifically to the β subunit of tubulin present on microtubules (Oberhoff et al., 2002;Axel et al., 1997); predominantly exerting its effect at the G 2 /M phase of the cell-cycle, but is also known to inhibit cell division at G 0 /G 1 (Oberhoff et al., 2002;Waugh and Wagstaff, 2004;Blagosklonny et al., 2004;Scheller et al., 2003). SIR, on the other hand, binds to cytosolic protein FKBP12, subsequently hindering the functionality of FKBP-rapamycin-associated protein (FRAP), interfering with the biological cascade occurring in the mTOR pathway as a result of growth factor stimulation. The rapamycin-FKBP12 (rFK) complex halts cell progression early, at the G 1 /S transition period through the synthesis of cyclin dependent kinase inhibitors, even in growth stimulated cells (Oberhoff et al., 2002;Wessely et al., 2006;Edinger et al., 2003).
Various studies have shown that these drugs inhibit cell growth in vitro. These include the demonstration of dose-dependent inhibition of cell growth by PTX for various exposure times (Axel et al., 1997;Scheller et al., 2003;Wiskirchen et al., 2004;Blagosklonny et al., 2004), where SMCs tend to be apoptosis-reluctant as PTX induces cell senescence, an irreversible arrest of the cell cycle through the stimulation of the p 53 tumour suppressor gene (Blagosklonny et al., 2004;Blagosklonny et al., 2006).
Similarly, literature illustrates that SIR (also known as rapamycin) markedly inhibits the growth and migration of SMCs, where the drug acts to prolong the G 1 phase of the cell cycle (Marx et al., 1995;Poon et al., 1996). Unlike PTX, cells inhibited by SIR retain their reproliferative potential at high doses (Blagosklonny, 2018). Further, Richard Gallo et al. (1999) were able to demonstrate that this inhibition occurred through increased levels of cyclin-dependent kinases. Lastly, Matter et al. (2007) compared SIR with tacrolimus (TAC), exhibiting that the former was more potent at inhibiting SMC growth and migration, while the latter was less potent at inhibiting endothelial cell (EC) growth, the fundamental issue in delayed arterial healing (Matter et al., 2007).
As improvements in ISR and ST rates are desired, other drugs have been studied besides PTX and SIR. Marra et al. (2000) considered the inhibitory effects of sodium salicylate (SAL) on cell proliferation, known to suppress cyclooxygenase-1 and 2 (COX-1, COX-2) transcription (Kenneth, 2000), where the latter is a significant mediator in the inflammatory response (Deng et al., 2001), inducing G1 cell-cycle arrest (Perugini et al., 2000). Recently, Sun et al. considered elemene, a drug which influences cell proliferation at the G 2 /M phase, similar to PTX (Sun et al., 2017). They presented positive results, where inhibition of SMC proliferation and migration was observed, whilst EC recovery was stimulated. Further analogues of SIR, such as everolimus, zotarolimus and biolimus were also shown to inhibit the growth and migration of SMCs in vitro (Lavigne et al., 2012;Chen et al., 2007;Kim et al., 2018).

Selection of data sets for mathematical model calibration and validation 2.2.1. Overview of in vitro data
Whilst there are several data sets in the literature that consider the effects of various drugs on SMC growth, these are not always presented in a format that enables mathematical model calibration and quantitative validation. For example, it is common for the action of drug on cell growth to be presented as a percentage of the control value, which is not stated (Axel et al., 1997;Matter et al., 2007;Lavigne et al., 2012;Chen et al., 2007;Zhu et al., 2009), or for information such as the initial and/ or final cell density to be omitted. Thus, the effect of the drug may only then be considered in a relative manner to the control case. Moreover, not all studies consider the effect of different drug doses, nor the influence of drug exposure time while these aspects are likely important in the context of arterial healing following stenting (Zhu et al., 2009;Sun et al., 2017;Kim et al., 2018). With this in mind, data sets that detail the effect of drug as % control and/or which do not consider the influence of multiple drug doses or drug exposure times have been neglected (Axel et al., 1997;Matter et al., 2007;Lavigne et al., 2012;Chen et al., 2007;Marx et al., 1995

Data selected for analysis
The data from two studies, those by Scheller et al. (Fig. 2) (Scheller et al., 2003) and Marra et al. (Fig. 3) (Marra et al., 2000) are considered for model calibration and validation. They consider two different drugs, PTX and SAL, which are reported to alter cell kinetics through intracellular receptors. These data sets satisfy our selection criteria such that: (i) data is presented as cell number/density, and not as '% control', absorbance, or % DNA synthesis with no reference to cell counts; (ii) multiple drug doses are considered; and (iii) drug retention is explored. However, we note that utilising data from two different studies does have some drawbacks, most notably in the interpretation of the results, where different experimental protocols and cell lines have been used. On the other hand, considering two different drugs holds the advantage of enabling us to consider whether the models presented are able to capture the data associated with drugs that have distinct physiochemical properties. Marra et al. (2000) studied the effects of SAL, a drug that has been previously utilised on cardiovascular devices (Jabara et al., 2008;Lee et al., 2014;Cysewska et al., 2019), and demonstrated the effect of three different doses of SAL on SMC proliferation for a 12 day incubation period. The authors also studied the effect of drug retention and subsequent cell recovery after removing the highest dose of drug following four days of exposure. Scheller et al. (2003) considered two doses of PTX within a neutral contrast media, exposed briefly to cells for 3, 10, and 60 min, where the supernatant was removed and SMC recovery was assessed across a 12 day period.

Modelling cell growth in the absence of drug
We start by providing a model for the control case, where cells are allowed to grow in the absence of drug, limited only by space constraints. We describe the number of cells as a function of time, c(t), via the logistic growth equation: where g is the intrinsic growth rate, K is the carrying capacity and c 0 the number of cells initially. Both g and K are taken to be constant: the former assumes that the cells are exposed to a constant growth stimulus for the duration of the experiment, while the latter may be related to the size of the culture plate.

Modelling the effects of drug exposure and retention on cell growth
We consider a series of models (D1, D2 and D3) of varying complexity to assess their ability to capture the effect of constant drug exposure on cell growth and, more importantly, the effect of drug retention on cell recovery (i.e. when the drug source is removed). Parameters that detail the action of drug (e.g. how quickly it (un) binds) are assumed to be drug dose independent.

Model D1
In our most sophisticated model, we assume that the level of drug bound to specific receptors dictates its effectiveness. Drug binding is modelled through a subtle amendment to the state-of-the-art nonlinear binding kinetics (Tzafriri et al., 2012), while (3.1) is modified to account for the effect of bound drug on cell growth: is the bound drug concentration as a function of time; k f and k r are the binding-on and binding-off rates, respectively; c d is the concentration of drug in the media, assumed to be homogeneously distributed and time-independent and; B is the average receptor density per cell. While the state-of-the-art nonlinear binding kinetics model (Tzafriri et al., 2012) assumes a constant density of binding sites through time (i.e. no proliferation of cells), the present model allows the maximum binding site density, Bc(t), to vary as a function of time in proportion to the current number of cells, i.e. when b(t) = Bc(t) all binding sites are occupied by drug, in other words, the receptors are saturated. Since receptor saturation has been linked to efficacy (Tzafriri et al., 2012), in (3.3) we assume the drug exerts its maximum effect on cell growth when bound drug levels are at their greatest, i.e. when b(t) = Bc(t) growth is completely halted. On the other hand, when b(t) = 0, the drug has no effect and growth is assumed to be logistic as in (3.1). All parameters are detailed in Table 1.

Model D2
Model D2 assumes instead that drug undergoes a simple cell internalisation process and elicits its effect through Hill kinetics, which have been used extensively in the literature to describe the dose-dependent effect of drugs (Gardner, 2000;Tham et al., 2008;Goutelle et al., 2009;Clairambault, 2009;Bernard et al., 2012;Ribba et al., 2014;Clarelli et al., 2020): In (3.4), b I (t) represents the internalised bound drug concentration as a function of time; α is the internalisation rate and; K P is the nondimensional partition coefficient such that, in equilibrium, b I = K P c d . In (3.5), k max is the maximum effect the drug has on cell growth and k 50 is the drug concentration at which the effect is half-maximal. This term establishes a non-linear relationship between b I (and the applied dose c d , via (3.4)) and its influence on cell growth, where the parameters k max and k 50 are dictated by the specific drug considered. We observe that when b I (t) = 0, the drug has no effect and growth is assumed to be logistic as in (3.1). All other parameters are as before, listed in Table 1.

Model D3
Our simplest model completely neglects drug internalisation and binding processes and replaces b I (t) in (3.5) with the applied dose, c d . Fig. 2. PTX at 14.6 μM exposed to bovine SMCs for 3, 10, and 60 min, and for 60 min at 1.46 μM. Acquired from (Scheller et al., 2003) with permission as it satisfies our selection criteria. This model represents a special case of (3.4-3.5), obtained by letting α→ ∞ and K P = 1. Whereas models D1 and D2, by their design, are able to account for the effects of drug retention after removal of the drug source, selecting c d = 0 in this model would simply return the logistic growth Eq. (3.1), meaning that the drug has no lasting effect. To account for this anomaly, we assume that when drug is removed from the media at time t = τ, the drug concentration associated with effect decays exponentially with associated decay rate k d . The model can then be separated into two parts, with (3.6a) describing cell growth while in the presence of drug and (3.6b) describing subsequent cell growth when drug has been removed:

Numerical solution of mathematical models
Prior to solving numerically, models D1, D2 and D3 were nondimensionalised for numerical convenience. For full details of the non-dimensionalisation procedure we refer the reader to the Supplementary Data (Section 1). All of the ordinary differential equations (ODEs) were solved using MATLAB's built-in 'ode15s' routine (Shampine and Reichelt, 1997), which integrates the equations over a specified time range with user-defined initial conditions.

Inverse estimation of model parameters
The data in Marra et al. (2000) comprises measurements of the total number of cells at various time points during no drug exposure (C m 12 ), or a constant exposure of either a low (L m 12 ), medium (M m 12 ) or high (H m 12 ) dose of drug. In an additional experiment (H m 4 ), the high dose of drug was removed after day 4 and subsequent cell growth was measured. These data are summarised in Table 2.
Unlike Marra et al. (2000), the experiments performed in Scheller et al. (2003) do not expose the cells to drug for a prolonged period of time. Instead, cells are exposed to a low dose of drug for 60 min (L s 60 ) or a high dose of drug for either 3 min (H s 3 ), 10 min (H s 10 ) or 60 min (H s 60 ), with subsequent cell proliferation assessed up to 12 days. In a control Fig. 3. A: Human saphenous SMCs exposed to 3 doses (0.1 mM, 1 mM, and 5 mM) of SAL for 10 days. B: 5 mM of SAL removed at day 4, illustrated by the arrowhead, depicting cell recovery upon the removal of drug. Experimental data reproduced from (Marra et al., 2000) as it satisfies our selection criteria.

Table 2
Summary of the five data sets provided by Marra et al. (2000). The first three data sets involve constant exposure of either a High (5 mM), Medium (1 mM) or Low (0.1 mM) dose of SAL for a period of 12 days. The fourth data set involves constant exposure of the high dose of SAL for 4 days, before the drug is removed from the medium, the cells washed, and the cells allowed to recover for the remaining 8 days. The final data set involves the control case where the cell growth in the absence of drug is measured over a period of 12 days.  Scheller et al. (2003) also consider cell proliferation in the absence of drug for a period of 12 days (C s 12 ). These data are summarised in Table 3. Scheller et al. (2003) do not provide data on number of cells and instead measure cell density. In the absence of any information in the paper regarding the surface area of the plates used, we assumed 48 well plates and converted the cell density data at each time point to cell number, for consistency with the Marra et al. (2000) data.
It is important to note that Scheller et al. (2003) do not provide any measurement data during their short drug exposure times (3, 10 or 60 min), with the first measurement taken after 3 days. This poses an additional challenge in terms of estimating the model parameters. Our approach is to run a forward simulation for the period of drug exposure, then export the number of cells and bound drug concentrations (for models D1 and D2) at the end of this period. These then form the new 'initial conditions' for the subsequent simulation of cell recovery in the absence of a drug source, where we then select c d = 0 in the models.
We implement a least squares approach within Matlab's 'fminsearch' routine and minimise the error (ε) given by the sum of the squares of the difference between model simulated data points (c j ) and the experimental data points (c e j ): where t i , i = 1,2,…,T, are the measurement time points. The subscript j, j = 1,2,…,N, denotes the number of data sets considered as part of the fitting process, with N the total number. When j = 1 we refer to an 'individual fit', whereas when j > 1 we refer to a 'multi-fit'. The 'fminsearch' routine adopts an error minimisation approach based upon the simplex method, such that the simplex consists of n +1 points, where n is the dimension of the problem. The routine executes various transformations (contraction, reflection, etc.) such that the final simplex comprises of parameters that present the smallest error (ε) within some tolerance limit (10 − 4 ). In the absence of reliable literature estimates of the model parameters, we adopt an unconstrained parameter estimation approach, i.e. we do not stipulate that the best-fitting parameters have to lie within a particular range.
Since the initial number of cells, c 0 , is known from the experimental data sets, the cell growth model (3.1) contains only two parameters, g and K, to be inversely estimated. For each study (Marra et al. (2000) and Scheller et al. (2003)) we inversely estimate the values of g and K such that the growth model given by (3.1) best fits the experimental data. The best-fitting parameters are then utilised in the subsequent models of cell growth subject to drug exposure (models D1, D2 and D3) which require varying numbers of additional parameters to be inversely estimated. In order to test the robustness of the models and parameter fitting procedure, we perform individual fits and multi-fits as well as predictions, as summarised schematically in Fig. 4.

Comparison between models and Marra et al. (2000) data
Our results demonstrate that the logistic growth model (3.1) is capable of capturing cell growth in the absence of drug (Fig. 5a). The best fitting intrinsic growth rate and carrying capacity were found to be Table 3 Summary of the five data sets provided by Scheller et al. (2003). The first three data sets involve exposure of a High (14.6μM) dose of PTX for 60 min, 10 min, and 3 min, respectively, followed by no drug exposure up to 12 days. The fourth data set involves constant exposure of the Low (1.46μM) dose of PTX for 60 min, followed by no drug exposure up to 12 days. The final data set involves the control case where the cell growth in the absence of drug is measured over a period of 12 days.  Fig. 4. Schematic of the parameter fitting process, broken into 3 main steps.
Step 1: Fit the growth model (Eq. (3.1)) to the control dataset, where parameters g and K are acquired and inputted as fixed parameters into step 2 and 3.
Step 2: Fit to each of the drug doses in each dataset separately (Eq. (3.7), j = 1) and subsequently predict the other responses based on the best-fitting parameter values obtained.
Step 3: Fit multiple data sets simultaneously (Eq. (3.7), j > 1) to attain best-fitting parameters, and use these to make predictions of cell growth for a different drug dose. This process is performed for each study (Marra et al. (2000) and Scheller et al. (2003)) separately. g = 0.204 s − 1 and K = 130, 777, respectively, with the error between the model and the data of ε = 0.0301.
We then proceeded to test the ability of the drug exposure and retention models (D1, D2 and D3) to capture the Marra et al. (2000) data when the cells were exposed to drug for the duration of the experiment. We observe in Fig. 5b-d that each of the models are able to capture the data well, when the parameter estimation is performed using each data set (H m 12 , M m 12 and L m 12 ) and model (D1, D2 and D3) separately. While the errors for each drug dose are found to be similar across the models, it is notable that using this method of fitting to the data sets separately, the best-fitting parameters vary not only across the various doses for a given model, but also between the models (Supplementary Data Table S1). This could suggest that the model parameters are in fact drug-dose dependent and not independent of drug dose as had been assumed in the model formulation. Other possible explanations could be that the parameters are not uniquely identifiable given the parameter estimation methodology we have used, or alternatively, that the models are unsatisfactory.
To probe this further, we then considered the case where the High drug dose is removed from the media after 4 days and the cells are allowed to recover (H m 4 ). For Models D1 and D2, we utilised the bestfitting parameters for each model based upon our fits to the H m 12 data set, and implemented c d = 0 from day 4 onwards to predict the cell recovery ( Fig. 6a and 6c). Our results demonstrate that the prediction using model D1 is extremely poor (ε = 0.6746), with the cell recovery severely underestimated. The prediction from model D2 is, however, substantially better (ε = 0.0603). One possible explanation for the poorer performance of model D1 is that the unbinding rate parameter (k r ) had been estimated under constant exposure conditions, where it may be expected that binding-on dominates. We therefore amended our parameter estimation approach by finding the model parameters that best-fit the following two data sets simultaneously: H m 12 and H m 4 . As can be seen in Fig. 6b this results in a substantially better fit (ε = 0.0125) for model D1 to data set H m 4 without adversely affecting the fit to data set H m 12 (ε = 0.0259). This approach was repeated for model D2 (Fig. 6d), resulting in a modest improvement in (the already good) agreement with the data (ε = 0.0161). Given the nature of model D3, predictions are not possible without the further estimation of the parameter associated with cell recovery (k d ). In this case, we utilised the best-fitting parameters based upon our fit for model D3 to the H m 12 data set and performed a further parameter estimation step to obtain the value of k d that best fits data set H m 4 , whilst all other fitted parameters (k max ,k 50 ) are kept fixed. Fig. 6e demonstrates that this model is able to capture the data very well (ε = 0.0358). Unfortunately, Marra et al. (2000) do not provide any further data where the drug source is removed after a period of time for other drug doses, to enable us to make any more predictions to further validate the models.
We next sought to test the ability of the models D1, D2 and D3 to predict cell growth under constant exposure of drug. Our aforementioned results suggest that the quality of the prediction will depend heavily on the data set(s) initially used to obtain the best-fitting parameters of the respective models. For each model, we firstly utilised the best-fitting model parameters for the L m 12 data set and then adjusted the drug dose parameter c d to predict the effect of medium and high drug doses on cell growth. This was then repeated for the case where the model parameters were obtained from the fit to the M m 12 data set (with the response to high and low drug predicted) and finally where the model parameters were obtained from the fit to the H m 12 data set (with the response to medium and low drug predicted). As expected, the quality of predictions were generally poor (Supplementary Data Table S2), owing to the wide variation in best-fitting parameters, depending on the data set used to perform the fit. In order to improve the predictions, we subsequently performed a series of multi-fits, where multiple data sets were used as part of the parameter fitting procedure, and then the model used to predict a remaining data set. We found that the best predictions using models D1 and D2 were obtained when the M m 12 , H m 12 and H m 4 data sets were used simultaneously to find the best fitting model parameters, with  Table S3). The best prediction using model D3 was obtained when the H m 12 and M m 12 data sets were used simultaneously to find the best fitting model parameters, with the effect of the low drug dose predicted (Fig. 7c) (Supplementary Data Table S3). As before, model D3 still required a subsequent fitting step using the H m 4 data set to estimate the parameter (k d ) associated with drug retention. As can be observed in Fig. 7a-c, model D1 performs poorest, while models D2 and D3 generally perform well.

Comparison between models and Scheller et al. (2003) data
Our results demonstrate that the logistic growth model (3.1) is capable of capturing cell growth in the absence of drug (Fig. 8a). The best fitting intrinsic growth rate and carrying capacity were found to be g = 0.5074 s − 1 and K = 86, 238, respectively, with the error between the model and the data of ε = 0.00382.
We then proceeded to test the ability of the drug exposure and retention models (D1, D2 and D3) to capture the Scheller et al. (2003) data (H s 60 , H s 10 , H s 3 and L s 60 ). We observe in Fig. 8b-d that each of the models are able to capture the data well, when the parameter estimation is performed using each data set (H s 60 , H s 10 , H s 3 and L s 60 ) and model (D1, D2 and D3) separately. Similarly to the Marra et al. (2000) data sets, we observe that while the errors associated with the fitting to each data set are found to be similar across the models, the best-fitting parameters vary depending on the data set and the model used (Supplementary Data  Table S4). This could again suggest that the model parameters are in fact drug-dose dependent; the parameters are not uniquely identifiable given the parameter estimation methodology we have used, or; that the models are unsatisfactory.
We then tested the ability of the models D1, D2 and D3 to predict cell growth. For each model, we firstly utilised the best-fitting model parameters for the H s 60 data set and then adjusted the exposure time to predict the shorter exposure cases (H s 10 and H s 3 ) before adjusting the drug dose parameter (c d ) to predict the L s 60 case. This was then repeated for each data set in turn, with the best-fitting parameters based on a given data set used in the predictions of the remaining cases. As expected, the quality of predictions were generally very poor (Supplementary Data Table S5), owing to the wide variation in best-fitting parameters, depending on the data set used to perform the fit. In order to improve the predictions, we subsequently performed a series of multifits, where multiple data sets were used as part of the parameter fitting procedure, and then the model used to predict the remaining data sets. In particular, we considered the following three cases, with the data sets used as part of the multi-fit specified: (i) H s 60 and L s 60 (Supplementary  Data Table S6); (ii) H s 60 and H s 3 and; (iii) H s 60 , L s 60 and H s 10 We found that the best predictions were obtained for case (i) and model D3 where the 60 min exposure of low and high drug dose data sets were used simultaneously to find the best fitting model parameters, with the effect of the shorter exposure times of the high drug dose predicted (Fig. 9c). As can be observed in Fig. 9a-b, even with this multi-fit approach, models D1 and D2 still result in a very poor prediction of the shorter high drug dose exposure times.

Fig. 7.
Comparison between data (dots) provided by Marra et al. (2000) and mathematical models (curves) of cell growth and recovery when a multi-fit is used to obtain the best fitting model parameters before the effect of the low drug dose is predicted (illustrated by the dashed line). Here, c 0 = 17, 100. (a) Model D1 (Eqs.

Discussion
The first two drugs coated on commerical DES, PTX and SIR, were already being used in other applications and were not designed specifically to counteract restenosis. Drug dosing for these early stents was based upon in vitro cell-based studies and the effectiveness observed in animal studies that followed. The subsequent success in humans at mitigating against restenosis has led to a general acceptance that the drug types and doses used on commercial DES are adequate, with the community predominantly focussing subsequent research efforts on new stent materials and geometrical designs. PTX, SIR and, more commonly, analogues of SIR are still the preferred drug choice for DES and there is only a modest variation in drug dose and release kinetics provided by existing commercially available stents. The current restenosis rates, while low, still equate to hundreds of thousands of patients worldwide, and with only modest improvement in patient outcomes despite the astronomical investment in innovative stent designs over the years, it is certainly worth revisiting the drug aspect of DES.
Surprisingly, little attention has been devoted in the literature to mathematically modelling the effect of drug on SMC proliferation, despite the fact that this process is a key driver of restenosis. A useful starting point is to develop mathematical models based upon in vitro data: while this represents a hugely simplified setting, any model which shows predictive power in this setting can be taken forward to the more complex in vivo situation with some confidence. To our surprise, while there exists several in vitro experimental studies in the literature, few provide sufficient details to enable quantitative comparison with mathematical models. In this paper, we identified two appropriate existing experimental studies and tested the ability of three mathematical models of varying complexity to describe their cell proliferation data. We now summarsie the key findings.

The logistic growth model well-describes cell proliferation in the absence of drug
Literature suggests that the logistic growth model is adequate in portraying cell growth in vitro (Jin et al., 2017). For each data set considered, Marra et al. (2000) and Scheller et al. (2003), we have demonstrated that the logistic growth model captures the cell proliferation data very well (Figs. 5a and 8a). While neither Marra et al. (2000) nor Scheller et al. (2003) provide any details related to the size of culture plates used, we note that our with inferred carrying capacities (K) are consistent with the range of cell capacity for most commonly used cell culture plates (Useful Numbers for Cell Culture thermo fisher, 2020). Our results illustrate a clear difference in cell growth between the two data sets. It is known that cell proliferation from species-to-species is  Table S4. ε H60 ε H10 , ε H3 , and ε L60 correspond to the error associated with the fit for each drug exposure data set ( different, which could provide justification to the fact that the growth rates (g) and K vary between the studies (Schwartz et al., 1996;Charlebois and Balázsi, 2019). Moreover, these differences may be explained by variations in experimental protocols (e.g. the size of culture plate used, not reported). While the Scheller et al. (2003) data appears to be reaching a clear plateau, the same cannot be said for the Marra et al. (2000) data: experimental data at further time points would shed light on this. If a plateau has indeed been reached, as Marra et al. (2000) suggest, then our inferred value for the carrying capacity may be an over-estimate.

Each model is capable of capturing all of the drug exposure data sets, but predictions are generally poor
We have demonstrated that when parameter fitting is performed on each data set and model separately, the agreement between each model and data set is generally very good (Fig. 5b-d and Fig. 8b-d), seemingly suggesting that each model is adequate. However, probing further, we found that the best-fitting model parameters varied widely across the different data sets (Supplementary Data Table S1 and S4). As a consequence, accurately predicting any given data set based on model parameters fitted to another data set, proved challenging. Indeed, we were only able to make acceptable predictions when multiple data sets were used as part of the fitting process. Such an approach may be justified by the fact that certain model parameters (for example unbinding rates) may not be uniquely identifiable given the available data. For example, in the case of constant exposure, binding-on may dominate over bindingoff. A key disadvantage of the multi-fit approach is that it then limits the number of available data sets available to predict. Our results suggest that our model parameters may well be drug-dose dependent and future experiments should be geared towards testing this hypothesis.

The simplest model out-performs the state-of-the-art non-linear binding model
For each data set, Marra et al. (2000) and Scheller et al. (2003), our simplest model D3 performs best in terms of both the goodness of fit to the data and also the ability to make accurate predictions, albeit after model parameters have been obtained following a multi-fit procedure (Figs. 7c and 9c). We now explain why this may be the case.
It is notable that the Scheller et al. (2003) data indicates no significant difference in cell proliferation regardless of whether the cells are exposed to the high dose of PTX for 3 min, 10 min or 60 min. In other words, even with a short exposure time to the high dose of PTX, the effect of the drug is long-lasting. To probe this further, we plotted the corresponding normalised bound drug concentrations for models D1 and D2 over the first 60 min (Fig. 10a-b). Fig. 10a-b demonstrate clearly that model predictions of normalised bound PTX concentration are different at the end of the different drug exposure periods for the high PTX dose (3 min, 10 min and 60 min). Since the effect of the drug on cell proliferation is directly related to the bound drug concentrations in model D1 and model D2, it therefore follows that these models are incapable of capturing the similar response to these different drug exposure times that is observed experimentally. More specifically, for Model D1 to be able to capture these data we would require a similar level of receptor saturation to be achieved after drug exposure times of 3 min, 10 min and 60 min. For this to be achievable and at the same time capture the effect of the low dose drug, we would require the model parameters to be drugdose dependent. While reliable estimates are unavailable in the literature for the majority of the parameters of our models, there are some estimates of k f and k r available for PTX (Tzafriri et al., 2009;Indolfi et al., 2016;Escuer et al., 2020). We have verified that using these values in model D1 results in a response indistinguishable from the control case, further questioning the applicability of the saturable receptor binding model to in vitro data. Model D3, on the other hand, completely neglects the drug binding process, with the effect of the drug instead related to the applied dose c d . This likely explains why model D3 outperforms models D1 and D2.
Turning to the Marra et al. (2000) data, we note that the low dose of SAL has a very small effect on cell proliferation, making it difficult to distinguish form the control (no drug) case, while the high dose of SAL almost completely inhibits cell growth. Estimating model parameters based on fitting to these data-sets and then making predictions of the remaining data set (the medium dose) is therefore less than ideal. It is notable that when SAL is removed from the media, the cells recover at a rate comparable to the control case, indicating that, unlike PTX, hydrophilic SAL has minimal lasting effect on the cells, suggesting that it not as strongly retained as PTX and that a dynamic binding/unbinding model may not be suitable. Fig. 7a shows that Model D1, even when a multi-fit is performed using the medium and high SAL data set, is unable to simultaneously capture the medium SAL dose data set and predict the low SAL dose data set. Models D2 and D3 neglect the drug binding process and any associated delay before the drug has any effect, which likely explains the better performance of these models for this data set.

Implications for future computational modelling and experiments to improve DES treatment
Our results put into question the applicability of the existing state-ofthe art model (D1) of drug kinetics following stenting, which assumes that efficacy is determined by nonlinear binding to cell receptors. While this model may be applicable to certain drugs, we have shown here that it is not capable of describing the effect of varying doses and exposure times of SAL or PTX on SMC proliferation in vitro. It is desirable to extend our analysis to sirolimus and analogues, which have become the drugs of choice for DES, but for the reasons mentioned in Section 2.2, we have been unable to identify any published studies involving these drugs where the data is presented in a format that readily allows comparison with mathematical models of the type presented here. Notwithstanding, PTX in particular is still widely used, especially on drug-coated balloons (DCB) for use in peripheral artery disease and coronary artery disease. With the safety of PTX recently being brought into question, particularly with respect to peripheral artery disease (Katsanos et al., 2018), the need to better understand the mode of action of this drug has never been more urgent. It may be that different models are required for different drugs, as has previously been suggested in the literature (Jarrett et al., 2019). It is also worth noting that our models do not include equations for nonspecific binding to tissue components, given the focus here is on in vitro studies. It may well be that drug effectiveness is dictated by some combination of specific and non-specific binding, which could at least partially explain any discrepancy with existing models.
While models of drug uptake and drug binding are typically parameterised using in vitro or ex vivo data, they are frequently utilised within more complex multiphysics models aimed at simulating the in vivo environment, where other effects such as blood flow, drug transport (from DES or DCB) and arterial healing are considered. Our results highlight that inversely estimating model parameters based on one dose of drug, then using the model to make predictions related to other drug doses without further validation may lead to inaccurate predictions. Thus, we strongly recommend parameterising models of the effect of drug on cell proliferation based on experimental data across multiple drug doses and exposure times, before incorporating within more complex models. In order to optimise the release of drug from DES, it is important to develop a clear understanding of the effect of different drug doses and exposure times on SMC proliferation since the in vivo situation dictates that SMCs will be exposed to spatiotemporal concentrations of drug, i.e. cells at different spatial locations will be exposed to different concentrations of drug at different times. Taken together, our results suggest that further experimental research, in tandem with mathematical modelling, is required to better understand the influence of varying doses of drug on cell proliferation, and ultimately the arterial healing response. Future experimental and modelling work should seek to address some of the limitations highlighted in this work by considering the possibility of model parameters being drug dose dependent; the cellcycle specific nature of the drugs and; introducing spatiotemporal gradients of drug. Ideally, details of future experimental protocols should be provided to a level that allows quantitative comparison with mathematical models, including information on the size of culture plates used as well as measurements on cell number or cell density provided. Finally, it is of interest to pursue reliable estimates of the various model parameters. This would enable a more robust parameter-fitting approach through the imposition of constraints and ultimately aid in uniquely identifying the various model parameters.

Conclusions
In this work we have devised mathematical models of varying complexity to probe their ability to capture in vitro SMC proliferation subject to drug exposure. Our results question current thinking about the mechanisms through which drugs coated on stents attenuate SMC proliferation, a key component of the arterial healing process. Future research should focus on assessing the translatability of cell based studies and associated mathematical models, to the reality of the in vivo situation. Such research should ideally involve mathematical modelling and experiments in tandem.

Declaration of Competing Interest
Keith Oldroyd is a full-time employee of Biosensors International.