Design Optimization of Tumor Vasculature-Bound Nanoparticles

Nanotherapy may constitute a promising approach to target tumors with anticancer drugs while minimizing systemic toxicity. Computational modeling can enable rapid evaluation of nanoparticle (NP) designs and numerical optimization. Here, an optimization study was performed using an existing tumor model to find NP size and ligand density that maximize tumoral NP accumulation while minimizing tumor size. Optimal NP avidity lies at lower bound of feasible values, suggesting reduced ligand density to prolong NP circulation. For the given set of tumor parameters, optimal NP diameters were 288 nm to maximize NP accumulation and 334 nm to minimize tumor diameter, leading to uniform NP distribution and adequate drug load. Results further show higher dependence of NP biodistribution on the NP design than on tumor morphological parameters. A parametric study with respect to drug potency was performed. The lower the potency of the drug, the bigger the difference is between the maximizer of NP accumulation and the minimizer of tumor size, indicating the existence of a specific drug potency that minimizes the differential between the two optimal solutions. This study shows the feasibility of applying optimization to NP designs to achieve efficacious cancer nanotherapy, and offers a first step towards a quantitative tool to support clinical decision making.

Nanotherapy may constitute a promising approach to target tumors with anticancer drugs while minimizing systemic toxicity. Computational modeling can enable rapid evaluation of nanoparticle (NP) designs and numerical optimization. Here, an optimization study was performed using an existing tumor model to find NP size and ligand density that maximize tumoral NP accumulation while minimizing tumor size. Optimal NP avidity lies at lower bound of feasible values, suggesting reduced ligand density to prolong NP circulation. For the given set of tumor parameters, optimal NP diameters were 288 nm to maximize NP accumulation and 334 nm to minimize tumor diameter, leading to uniform NP distribution and adequate drug load. Results further show higher dependence of NP biodistribution on the NP design than on tumor morphological parameters. A parametric study with respect to drug potency was performed. The lower the potency of the drug, the bigger the difference is between the maximizer of NP accumulation and the minimizer of tumor size, indicating the existence of a specific drug potency that minimizes the differential between the two optimal solutions. This study shows the feasibility of applying optimization to NP designs to achieve efficacious cancer nanotherapy, and offers a first step towards a quantitative tool to support clinical decision making.
Targeted cancer nanotherapy relies on nanocarriers to deliver anticancer agents safely to tumors while minimizing systemic toxicity. Nanocarrier-mediated drug delivery has been associated with up to an 8-fold increase in drug efficacy compared to conventional chemotherapy 1 . Both drug and nanocarrier design play important roles in treatment efficacy. In general, drug design focuses on finding compounds that inhibit cancerous cell viability or proliferation, while nanocarrier design aims at developing nano-vehicle structures that maximize drug concentration in tumor relative to healthy tissue, thus reducing adverse drug effects. Experimental and computational methods have been employed to pursue such designs.
Specifically for vasculature-bound nanoparticles, in vitro studies have focused on characterizing the effect of nanocarrier design on margination, adhesion, and uptake while flowing through the tumoral vasculature. The tendency for vascular-borne liposomes and metal nanoparticles to drift from the blood streamlines towards the tumor vessel walls was studied by Toy et al. 2 using an in vitro microcirculation model. Considering different designs that vary in nanocarrier diameter {56, 60, 65, 100, 300 nm} as well as aspect ratio {0.45, 1}, it was shown that small eccentric nanoparticles are associated with stronger margination tendency. Larger nanoparticles, however, have exhibited a different correlation. Charoenphol et al. 3 examined the margination of nanoparticles of larger diameters {200, 500, 2000, 5000 nm}, flowing in an in vitro parallel flow chamber, showing that the margination rate increases with size. Patil et al. 4 measured in vitro the adherence strength of nanoparticles with diameters {5, 10, 15, 20 μm} coated with P-selectin glycoprotein ligands, showing that large nanoparticles have strong adherence properties due to a high contact area with the vascular endothelium. Importantly, large nanoparticles are subjected to stronger hemodynamic forces and torques that may dissociate them from vessel walls 5 . Boso et al. 6 measured the accumulation of nanoparticles of {0.75, 1, 2, 4, 6 μm} diameters in a parallel flow chamber. Data was fitted using an artificial neural network to correlate nanoparticle accumulation and size. As the size increases, the accumulation increases until it saturates or starts declining after a certain diameter that ranges between 4 and 6 μm depending on the wall shear rate. The study indicates the presence of moderate nanoparticle size that maximizes the adherence properties.
In addition to the aforementioned in vitro studies, in vivo investigations have evaluated the overall efficacy of nanoparticles in living subjects. Rostami et al. 7 showed that encapsulating Doxorubicin (DOX) in H6-equipped SCiEntiFiC REPORTs | (2018) 8:17768 | DOI: 10.1038/s41598-018-35675-y nanocarriers trebles the inhibition of mammary gland tumors in a mouse model compared to free DOX. Docetaxel-loaded nanoparticles of 349 nm diameter were delivered to mouse mammary tumors in 8 . Significant improvement in antitumor activity was obtained by delivering the drug through nanoparticles. The effect of nanoparticle size was studied by Joshi et al. 9 . Liposomes of diameters {60, 80, 200, 650, 670 nm} delivered to gliomas have indicated that 200 nm had the highest uptake rate. Other in vivo studies are reviewed in Zhang et al. 10 , for which nanoparticle design recommendations were based on increasing circulation time, taking advantage of the enhanced penetration and retention (EPR) effect, and maintaining high drug entrapment efficiency in the nanoparticle synthesis stage. Quantifying nanotherapy efficacy using in vitro assays may require lengthy preparatory steps, which include setting up proper cell lines and reagents, synthesizing nanoparticles, and tailoring experimental protocols. The complexity of these studies is further escalated in vivo. The cost and time associated with in vivo studies present limitations to evaluating different designs. Not only is acquiring and maintaining animal models expensive, but there may exist a long process from the initiation of oncogenic mutation or transplantation of xenografts, tumor proliferation, to monitoring tumor regression after nanoparticle injection. This also requires advanced imaging techniques and multidisciplinary expertise. For these reasons, computational modeling offers an attractive option for exploratory evaluation of nanoparticle design that complements experimental work, including investigation of a wide range of variables.
Computational modeling of tumor growth and nanoparticle delivery is, however, non-trivial. Such models need to consider a variety of biological processes such as angiogenesis and drug cellular uptake in order to yield informative results. The models typically consist of submodels, coupled sequentially or iteratively, which may be difficult to solve and slow to converge to a solution. Decuzzi et al. 11 modeled the margination of nanoparticles by taking into consideration buoyant and hemodynamic forces, as well as van der Waals interactions. The model explored a wide range of sizes, with results showing that diameters between 100 and 400 nm have relatively slower margination rates. Decuzzi and Ferrari 5 modeled the probability of nanoparticle adhesion at the vasculature wall. Parameters included nanoparticle size, aspect ratio, vessel wall shear stress, receptor-ligands association constant, ligand density, and receptor density. In parallel, tumor growth in two spatial dimensions coupled with neovasculature development was modeled mathematically in [12][13][14][15] . In 16 , Frieboes et al. integrated the nanoparticle delivery model of 5 with the tumor mechanics models of 14,15 to create a comprehensive model to predict the intra-tumoral distribution and accumulation of vasculature-bound nanoparticles. However, the computational cost of the integrated model hinders the evaluation of all possible nanoparticle designs of interest. Nanospheres of {100, 600, 1000} nm diameters were simulated using different values of nanoparticle avidity and tumor conditions. The results showed that large nanoparticles accumulate at higher rates at the tumor periphery, while smaller nanoparticles have lower adherence strength but are distributed more uniformly throughout the tumor tissue. Van de Ven et al. 17 used the model in 15 to study the effect of drug potency on tumor growth inhibition and to determine the number of nanoparticles needed to reach a half maximal effective concentration IC 50 . Wu et al. 18 used the model in 15 to study how the tumor interstitial pressure and fluid flow affect nanoparticle transport and distribution. Recently, the model in 16 was extended to simulate the tumor response to drug release from vasculature-bound nanoparticles 19 . Further information regarding mathematical characterizations of tumor nanotherapy can be found in 20,21 .
Previous computational and experimental work has primarily investigated the performance of nanocarriers using selected values of design variables such as nanoparticle size, aspect ratio, and ligand density. Optimizing nanoparticle design, however, requires design space exploration, which may be impractical to accomplish solely via empirical methods or computational models. Recently, Chamseddine and Kokkolaras 22 addressed this issue by applying rigorous optimization to the design of nanoparticles in order to maximize the tumoral nanoparticle accumulation and distribution with respect to the nanoparticle physical and chemical properties. The model was static, i.e., considered a single injection and does not update the tumor size and vasculature in response to the treatment. Although a sensitivity analysis was conducted to prove the robustness of the optimal design with 20% change in the average wall shear stress, the proposed design is not guaranteed to remain optimal if the tumor structure changes drastically. In this paper, optimization is applied to a "blackbox" version of the model, presented in 19 , that considers the dynamic changes in the tumor size and vasculature, enabling to obtain a nanoparticle design that is potentially optimal over the course of treatment. This represents a first step toward the goal of developing a clinically-relevant numerical tool to assist nanoparticle design on a patient-specific basis.

Computational Model.
A previously developed numerical model is used to compute tumor growth and response to nanoparticle drug delivery. The model 19 builds upon the work of 5,[13][14][15][16] , and is used in this study as a "blackbox" system with a limited set of inputs to calculate nanoparticle accumulation and tumor regression as a function of the nanoparticle design. Briefly, the model is composed of 4 submodels that are coupled in the configuration shown in Fig. 1a 16 . The "Tumor Compartment" submodel computes the progression/regression of the tumor as a function of drug, oxygen, and nutrient concentrations. It also creates hypoxic and necrotic regions, which induce angiogenic factors (TAF). TAF drives the "Angiogenesis" submodel to develop new blood vessels. "Angiogenesis" is coupled with the "Flow" submodel, which determines the wall shear stress among other flow properties in the preexisting and neo-vasculatures. The "Nanoparticle" submodel determines the nanoparticle accumulation and drug release to the cancerous tissue as a function of the nanoparticle design and tumor parameters. The model does not currently include an immune system response. Thus, the nanoparticles are essentially considered hidden from immunosurveillance. Table 1 lists the input parameters to the "blackbox" system and their values as used in this study, unless it is otherwise specified in the text. The complete list of the computational model parameters is reported in [14][15][16]19 .   Growth Phase. At day 0, a transformed group of cells is placed in the middle of a two-dimensional panel representing a tissue with blood vessels that are laid orthogonally as shown in Fig. 1b to simulate the regular vascularization of normal tissue. Blood enters the tissue from the bottom left corner to supply oxygen and nutrients. This enables the cancerous cells to proliferate and develop into a mass (tumor) due to the suppression of apoptosis. As the tumor grows, some cancerous cells distal from the blood vessels become hypoxic. This tissue produces TAF to stimulate new blood vessels that sprout from the existing vessels in order to supply the tumor tissue with blood. The neovasculature has an irregular structure and promotes tumor progression as shown in Fig. 1b. The growth phase is stopped at day 18 when the tumor reaches a diameter of 780 μm, after which the treatment phase starts.
Treatment Phase. Drug-carrying nanoparticles are injected into the blood vessel inlets at day 18. A fraction of the nanoparticles adhere to the tumor vessels. This fraction depends on the nanoparticle design. Anticancer agents are then released to the cancerous tissue. If the drug concentration exceeds a specific threshold, the tissue will die via apoptosis depending on the drug potency (see Table 1). A preliminary investigation of a 6-day treatment phase was simulated using different nanoparticle diameters. The change in tumor diameter in response is depicted in Fig. 1c. Since the curves do not intersect after 1 day of treatment, the treatment duration does not affect the relative performance for different nanoparticle sizes. Therefore, in our search for the optimal diameter, the treatment phase is stopped after 36 hours of nanoparticle injection, saving substantial computational time. Note that after a certain time of nanoparticle injection, a relapse is observed. This regrowth could be caused by two factors: either the drug is exhausted or the drug does not reach cytotoxic concentration for all of the proliferating tumor tissue.
Optimization. The ultimate goal of drug-based cancer treatments is to eradicate tumors completely or reduce their size prior to radical treatment intervention as a neoadjuvant therapy. Additionally, the aim of using nanoparticles as drug carriers is to reduce the side-effects associated with conventional chemotherapy while maximizing the drug delivery to the tumor tissue. Hence, we consider two objectives: minimizing the tumor size at the end of the treatment phase, and minimizing the treatment toxicity by maximizing the accumulation of nanoparticles in the tumor. Accordingly, two objective functions are defined: Let x denote the set of nanoparticle design variables. In general, the vector x may include nanopaticle diameter, aspect ratio, elasticity, ligand density, ligand-receptor affinity constant, drug release rate, and drug load. In this study, we consider spherical nanoparticles (aspect ratio of 1) because they are easy to manufacture 23 , and because they are predominant in current clinical and experimental studies 24,25 . For instance, the clinically proven nanodrug Doxil -used to treat different types of cancers such as breast and ovarian -is composed of 100 nm spheres. The optimal values of x that minimize TD, called minimizers of TD, can be obtained by solving: where n is the number of variables considered, and l b and u b are the lower and upper bounds that define the feasible region of x. Constraints are implemented implicitly in the analysis model. The maximizers of TNP are determined by solving: Problems (1) and (2) are solved using the Mesh Adaptive Direct Search algorithm (MADS) 26 . Since the gradients of the computational model cannot be approximated reliably, we use MADS, which is a derivative-free optimization algorithm that has rigorous convergence properties. Moreover, the computational model used here as a blackbox for analysis is computationally expensive. To examine a single nanoparticle design, an Intel(R) Core(TM) i7-3770 CPU @ 3.4 GHz processor requires 1.5 hours of CPU time. For this reason, we use a surrogate-assisted optimization approach to reduce the number of computational model evaluations required to obtain the optimal design. Specifically, we utilize the search step of the iterative mesh adaptive direct search (MADS) algorithm 27 to solve a surrogate of the optimization problem, i.e., we solve the optimization problem using surrogate models of the objective functions. To enhance the efficiency of the MADS algorithm we build and update ensembles of surrogates and use a novel order-based error metric tailored specifically for surrogate optimization to utilize the best surrogate at each iteration 28 . In this manner, we generate several candidates for the next iterate, which we combine with the candidates generated at the poll step of MADS, which is the foundation of its convergence properties. The computational model is then used only to evaluate all these candidates opportunistically to select the next iterate. In other words, we generate a lot of useful information by means of computationally inexpensive surrogate models but make algorithmic decisions using the high-fidelity computational model.

Results
Optimizing Nanoparticle Diameter. The drug biodistribution depends on the diameter d of the drug-carrying nanoparticles, which are localized depending on their size 16,22,29 . Small nanoparticles are associated with longer circulation time, increasing their chance to reach the tumor. Their adherence properties are poor, however, due to the low contact area between ligands on the surface of the nanoparticles and receptors over-expressed in the vascular endothelium of the malignant lesion. On the other hand, large nanoparticles have strong binding affinity, but they are also exposed to high hemodynamic loadings that may dissociate them from the endothelium. In addition, large nanoparticles tend to accumulate at the periphery of the tumor or bind to healthy tissues before they reach the tumor site due to their low circulation time. The optimal nanoparticle diameter d * lies within this range and is obtained by solving Problems (1) and (2) while setting x = [d]. Nanoparticles smaller than 10 nm are exposed to renal clearance, and nanoparticles larger than 1000 nm may not be able to flow in narrow tumor vessels and cause embolism. Therefore, d lies between l b = 10 nm and u b = 1000 nm in this study. For this evaluation, the nanoparticle avidity α = 2.95 × 10 10 m −2 , calculated by considering typical values for the receptor density, ligand density, and receptor-ligand binding constant under zero load, which are in the order of 10 12 #/m 2 , 10 14 #/m 2 , and 10 −14 m 2 respectively 16 . Figure 2a compares empirically-selected points with MADS-selected points in an attempt to minimize TD. In empirical methods, trial points are randomly chosen; however, using MADS, trial points are selected systematically to converge to the optimal solution. Note that when MADS approaches the optimal diameter, it tries many points in the vicinity before terminating at the best solution. The obtained optimal nanoparticle diameter up to 1 nm accuracy is d * = 190 nm, reducing TD * = 0.683.
Similarly, problem (2) is solved to maximize the tumor nanoparticle accumulation. The progress of MADS is shown in Fig. 2b. The optimal nanoparticle size that maximizes TNP is d * = 147 nm leading to TNP * = 0.137; i.e., around 14% of the injected nanoparticles successfully reach and adhere to the tumor. The corresponding TD is 0.484. Both optimal solutions are summarized in Table 2.  Effect of Nanoparticle Avidity. Vasculature-bound nanoparticles are equipped with ligands of high binding affinity to receptors over expressed in the vascular endothelium of tumor vessels. In the model, we assume that the integrin α ν β 3 exists with an area density m r in the malignant lesion. Corresponding ligands such as vitronectin 30 , fibronectin 31 , fibrinogen 32 , and osteopontin 33 are available at the surface of the nanoparticles with density m l . Each receptor-ligand pair has a certain affinity that is quantified by the binding constant under zero load K A 0 . The nanoparticle avidity α ∝  m m K r A 0 corresponds to the overall affinity of the nanoparticle 16 . The optimization problems (1) and (2) are solved again using a different value of α to check if the optimal diameters change. Changing the parameter α has an effect on the optimal nanoparticle diameters as shown in Table 3, where the cases of α = 1 × 10 11 m −2 and α = 1 × 10 12 m −2 are listed. It can be observed that as α increases, the diameter that maximizes TNP decreases. This decrease in d * can be explained by the increased ligand and receptor densities to the point that the nanoparticle-endothelium contact area required to cause binding is reduced. The minimal TD is altered, however, since smaller nanoparticles have lower drug loads. Since the value α has an effect on the optimal nanoparticle diameter, it is necessary to add it to the set of variables to optimize it along with d in an all-in-one optimization problem.
Treating Nanoparticle Avidity as a Design Variable. Let x = [d, α] T . The range of α is assumed to be between 10 10 and 10 12 m −2 complying with typical ranges of m r , m l , and K A 05, 16 . Hence, the feasible design space is l b = [10 nm, 10 10 m −2 ] T and u b = [1000 nm, 10 12 m −2 ] T . Results for minimizing the tumor diameter and maximizing the tumor nanoparticle accumulation are shown in Table 4.
Optimizing both nanoparticle diameter and avidity provides a better tumor reduction. Figure 2f displays the tumor at the beginning and end of the treatment, showing that the tumor reduces to 50.5% of its diameter at the start of treatment. Note that the optimal value of α lies at its lower bound (10 10 m −2 ) for both problems. The corresponding optimal diameter is increased to maintain an adequate contact area with the endothelium.
Relaxing the α-Boundary Constraint. The lower bound of α is an active bound; i.e., if it changes, the optimal value of α changes. To examine the proposed nanoparticle design practice of reducing the value of α and selecting a proper diameter, we reduce the lower bound of α to 10 9 m −2 and check if α * remains a boundary optimum. Figure 3 plots MADS progression toward the optimal solution of the relaxed problem. The solution confirms the existence of α * at the active bound. The variable d * remains an interior optimum having a value of 980 nm, increased to compensate for low ligand density.

Robustness of the Optimal Design.
The solution for problems (1) and (2) may change if the tumor morphology changes. Although numerical optimization can be a powerful tool that supports precision medicine dealing with patient-specific situations, designs that are aimed at treating a wide range of patients need to be insensitive to changes in tumor parameters. A rigorous method to attain robust designs is to optimize under uncertainty 34 , which will be addressed in future work. Here, we perform a sensitivity analysis with respect to model parameters that characterize the tumor microenvironment. We identify β and γ as candidates for altering the optimal design. The parameter β ∝ χμ/(k B Tm r ) combines Boltzmann thermal energy, blood viscosity, ligand-receptor binding force, and receptor density. The parameter γ is inversely proportional to the receptor density. More details about these parameters can be found in 5,16 .  (1) and (2) Table 3. Solutions of Problems (1) and (2) (1) and (2) The parameters β and γ are expected to have a direct impact on nanoparticle accumulation because they model nanoparticle-to-endothelium interactions. Therefore, we investigate the change in TNP with respect to the input vector using the interaction plot of Fig. 4. The interaction plot is a matrix plot, where the diagonal of the plot displays the categorical variables. The interaction of the parameter highlighted at the row-diagonal (i, i) with the parameter at the column-diagonal (j, j) is displayed at the off-diagonal position (i, j). For instance, the subplot (1, 2) plots the interaction of d and α. The horizontal axis is Z α , the vertical axis is the output TNP, and the different lines are the different values of Z d (indicated on the legend to the right of the corresponding row). The contribution of β and γ is illustrated in the subplots (1, 3), (1, 4), (2, 3), (2, 4), (3, 1), (3,2), (3,4), (4, 1), (4,2), and (4,3). In all of these plots, the graphs are either horizontal or coincide. Therefore, nanoparticle accumulation highly depends on the nanoparticle design and less on tumor biological conditions such as receptor density and blood properties.
Effect of Drug Potency. The potency of the drug encapsulated in the nanoparticles has an important role in tumor regression 17 . Drugs with higher potency cause fast shrinkage but are associated with high systemic toxicity. On the other hand, drugs with lower potency may evince slower tumor regression but have higher median toxic dose, lowering the associated adverse events. The computational model accounts for the drug potency through the proliferative term λ

ffect D D T drug
, which quantifies the interplay between cell mitosis, promoted by the availability of nutrients and oxygen σ, and cell apoptosis, which occurs if the drug concentration D exceeds a specific threshold T drug in the tissue 17,19 . The drug potency is measured by λ effect having a unit of effect per drug concentration. Clinically, drugs with higher potency have lower half maximal inhibitory concentration (IC 50 ). The parameter C D is a rescaling factor and A is the natural apoptosis rate.
In the Effect of Nanoparticle Avidity section, a hypothetical drug of moderate potency was used. The parameter λ effect was normalized to 1 for the drug considered. Drugs with higher potency are characterized by λ effect > 1, while drugs have a value of λ effect between 0 and 1. Although drug potency does not affect nanoparticle accumulation, it has an impact on the amount of drug needed to induce cell apoptosis, which is expected to change the minimizers of TD.
We showed before that there is a difference between the optimal diameter ⁎ d TD that minimizes TD and ⁎ d TNP that maximizes TNP. Figure 5a plots both optimal diameters for the case of α = 1 × 10 10 m −2 . If nanoparticles are smaller than =  Parametric Study with respect to Drug Potency. Figure 5 implies that the maximal therapeutic potential is not necessarily tied to the maximal nanoparticle accumulation. In fact, depending on drug potency, a smaller number of large nanoparticles can cause better tumor reduction than smaller nanoparticles. This defines the shaded zone of Fig. 5a. The case of a drug with lower potency (λ effect = 0.5) is then considered to investigate the change in the ⁎ d TD . For instance, assume a drug with λ effect = 1, used in previous sections, refers to Doxorubicin. If we consider the case of targeting invasive breast carcinoma, for which mutated MDA-MB-231 cell lines are a representative example with Doxorubicin IC 50 of 1.24 μM 35 , then λ effect = 0.5 refers to a drug of IC 50 = 2.48 μM, similar to Ponatinib 35 . The shaded zone becomes wider since ⁎ d TD is increased to 384 nm to provide higher drug volume needed in the tissue (Fig. 5b). In what follows, we assume λ effect = 1 refers to Doxorubicin and we use it as a benchmark to make all other comparisons with drugs of different values of λ effect .
In addition, a drug with higher potency is studied by setting λ effect to 2, simulating Lestaurtinib for example. Figure 5c shows that the width of the shaded region decreases. If λ effect increases further to 5, for instance if Midostaurin is used, the minimizer of TD becomes less than the maximizer of TNP. The reason is that a large drug load per nanoparticle is not needed to cause apoptosis at high values of λ effect . Therefore, the optimal solution shifts to small nanoparticles because they distribute more uniformly.
Notably, in all the considered cases, a small sacrifice in TNP leads to an increase in TD. Therefore, from a computational point of view, nanoparticle designs should be driven by minimizing TD. However, this conclusion may not be generalized; it requires extensive experimental support and should be evaluated for specific tumors. Furthermore, the reason TNP is less sensitive to the nanoparticle design could be due to an implicitly specified model parameter. Determining the drug potency that minimizes the discrepancy between the optimizers of the two objective functions would reduce the tradeoff between treatment speed and toxicity.
Conjugating ⁎ d TD and ⁎ d TD . In order to find a single nanoparticle diameter that optimizes both objective functions, we define the optimization problem Problem (3) aims at determining the drug potency at which the minimizer of TD coincides with the maximizer of TNP (288 nm). Solving problem (3) requires two loops. The inner loop computes ⁎ d TD given a drug potency λ effect that is specified by the outer loop. The outer loop iterates to minimize the difference between ⁎ d TD and 288 nm with respect to λ effect . In each outer loop iteration, the inner loop has to complete a full optimization process to find the minimizer of TD. The nested nature of problem (3) requires extensive computational time, which could exceed a month if the computational model used as a blackbox for analysis is used. Alternatively, a surrogate model is created by fitting the points that were evaluated earlier. Figure 6a shows the kriging metamodel constructed using the DACE (Design and Analysis of Computer Experiments) 36 . Exponential correlation functions and second-order polynomial regression models are employed to generate the kriging metamodel.
The solution process of problem (3) is illustrated in Fig. 6b. Each inner loop has a fixed value of λ effect on the vertical axis. Given λ effect , MADS visits the surrogate model and finds d that is closest to 288 nm, marked by the vertical line in the plot. The inner loop iterates horizontally to converge to the optimal solutions, shown in crosses. Then λ effect changes in the outer loop and the same procedure repeats until the optimal solution λ ⁎ effect is obtained. The minimal difference between ⁎ d TD and ⁎ d TNP is 6 nm. It is attained at the optimal drug potency λ = .

Multiple Treatment Cycles.
We consider the case where λ effect = 1 and α = 10 10 m −2 . A treatment with d = 334 nm showed a 50% reduction in tumor diameter, illustrated in Fig. 2f. In this section, we simulate three injections to reduce TD further assuming daily administration of nanoparticles. We obtained the optimal nanoparticle diameter for the second cycle of treatment by solving the optimization problem of equation (1)  , TD is further reduced to 0.22. Becoming very small, the treatment may be maintained in the third cycle using multiple nanoparticle diameters and not merely through the minimizer of TD. For this reason, the maximizer of TNP is chosen for the third injection to economize on the toxicity while achieving the same objective of continuous reduction in TD. The solution of TNP maximization in the third cycle is d* = 729 nm. The change in TD in the three treatment cycles is shown in Fig. 7.

Discussion
This study applies optimization to the design of drug-carrying nanoparticles targeting tumor vascular-endothelium. Empirical methods to obtain optimal designs are typically based on trial and error schemes. Instead, an optimization approach systematically converges to the optimal design using a significantly lower number of trial points. The reduction in computational time is necessary since optimal solutions are subject to change due to variations in tumor morphology, patient status, and cancer type 37 . Our previous study 22 was a first attempt to engage rigorous numerical optimization with the design of nanoparticles. We based our analysis on a model that simulates the distribution of nanoparticles and drug release to the tissue without considering the tumor response. In this study, we combined the optimization technique of 22 with a computational model of tumor nanotherapy 19 that models tumor evolution and includes complex biological processes such as angiogenesis and tissue necrosis [13][14][15][16] . This enables more reliable nanoparticle designs that are optimal for the whole duration of the treatment. First, the nanoparticle diameter was optimized. Previous studies have predicted that the optimal design moderately lies between lower, more uniform nanoparticle distribution, and higher less uniform distribution 16,19,38 . Experimental studies have examined a finite set of diameters to approximate the ideal size using in vitro and in vivo studies as reviewed in the Introduction and in 39,40 . Computationally, various diameters were investigated to correlate nanoparticle size and pharmacokinetics 29 . In 16 , three diameters were studied ({100, 600, 1000} nm). The respective tumor accumulation was {12.2%, 2%, 0.8%}. In comparison, the optimal diameter of 147 nm obtained in this study yielded accumulation of 13.7%.
The set of design variables has been expanded to include nanoparticle avidity because it has an impact on nanoparticle distribution 41 . Multidimensional search schemes give rigorous optimization a computational advantage over brute force methods that change one variable at a time. Optimizing nanoparticle diameter and avidity in an all-in-one problem shows a 1.6-fold decrease in tumor diameter and the same percent of nanoparticle accumulation as compared to optimizing the nanoparticle diameter alone. The resulting optimal diameters that maximize tumor targeting and minimize tumor diameter are 288 nm and 344 nm respectively. These diameters lie in the range that benefit from the EPR effect 6 .
The values of the optimal solution have a significant physical meaning. Previous studies have discussed that one way to prolong nanoparticle circulation is size reduction [42][43][44] . The results suggest that decreasing the nanoparticle avidity is a substitute approach in order to leverage the therapeutic potential of larger nanoparticles (of higher drug load). To further investigate, we reduced the lower bound of the nanoparticle avidity. The results confirm the existence of optimal nanoparticle avidity at the lower bound, whereas the optimal nanoparticle diameter increased to compensate for the decrease in the ligand-receptor pairing per unit area. This conclusion is specific for vasculature-bound nanoparticles where the drug is released at the endothelial layer. For the case where nanoparticles internalize to the tissue, smaller nanoparticles show higher cellular uptake rate, as reported in [45][46][47] .
In addition, small nanoparticles incur manufacturability limitations in terms of drug load since the efficiency of encapsulating free drugs depends on the nanoparticle size 8,[48][49][50][51][52][53][54][55] . For instance, we refer to the preparation of Doxorubicin-loaded albumin nanoparticles in 48 , as Doxorubicin is one of the most widely used antineoplastic agents 56 . The minimum nanoparticle size produced in 48 is 128 nm and was associated with 58% entrapment efficiency. However, based on our recommendation of using lower ligand density and moderate diameter of 334 nm, the entrapment efficiency could reach up to 78%. Hence, avoiding the unnecessary decrease in nanoparticle size may remove a fundamental constraint that hinders the development of efficient drug-loaded nanocarriers.
The robustness of the optimal design with respect to key tumor properties was then examined by studying the interaction among different parameters (Fig. 4). The results show that nanoparticle accumulation predominantly depends on the nanoparticle design rather than on tumor vessels properties such as receptor density, blood viscosity, and temperature.
The results indicate that maximizing nanoparticle accumulation is not optimal with respect to the tumor diameter reduction. Larger nanoparticles are needed to yield better tumor regression due to the associated high drug load. This was confirmed by running a study that uses a drug of lower potency. Results showed that the minimizer of tumor diameter is larger. Therefore, the optimal design is not unique; however, a range of optimal values can be recommended. This range is bounded by the maximizer of nanoparticle accumulation and minimizers of tumor diameter, where a tradeoff between the two objective functions exist. Depending on the patient status and clinical evaluation, an optimal solution can be selected to satisfy the weight given for each treatment attribution (fast tumor size reduction versus low systemic toxicity). In contrast, when the drug potency is high, less drug load per nanoparticle is required, and the width of the tradeoff range decreases. We increased the drug potency until the minimizer of tumor diameter is smaller than the maximizer of tumor accumulation. If the drug potency is high, smaller nanoparticles are expected to yield better tumor shrinkage because they distribute more evenly. These results are consistent with 16,19 .
The results motivated another study to find the drug potency at which minimizing tumor diameter is tied to maximizing nanoparticle accumulation. If this drug potency exists, nanoparticles could be designed to eliminate anticipated tradeoffs between the desired treatment attributes. Therefore, we formulated an optimization problem to find the drug potency that minimizes the difference between the two optimal designs. The complexity of the optimization problem leads us to synthesize a surrogate model using kriging metamodeling. The solution of the optimization problem indicates that a drug potency that is 4.6 times higher than the original value yields a minimal difference between the two optimal diameter values. A drug of this property would be similar to Omipalisib of IC 50 = 0.27 μM, shown to eliminate MDA-MB-231 cell lines 35 .
Integrating drug potency in the design of nanoparticles enhances drug targeting. Consider for instance the case of invasive breast carcinoma. Choice of targeted drugs such as small molecule inhibitors must be made depending on the driver oncogene responsible for the patient tumor. Talazoparib and Olaparib are two drugs used to target PARP1 and PARP2 genes that could be mutated in advanced breast tumors. Their IC 50 35 . This corresponds to λ effect of 3.9 for Talazoparib and 13.4 for Olaparib. Potentially, the IC 50 could be evaluated in the laboratory with the patient's tumor cells post biopsy, e.g., via monolayer or 3D cell culture cytotoxicity studies in order to help determine which drug is stonger. Toxicologically, Talazoparib represents the better option to target PARP1 and PARP2. If, however, the nanoparticle design that minimizes tumor diameter is predicted by the model to lead to low nanoparticle accumulation, i.e., wide grey region of Fig. 5, then minimizing TD with nanoparticles loaded with Olaparib may offer a better choice with nanoparticle diameters that are closer to the maximizer of TNP. We further used the optimization approach repetitively to explore the optimization of multiple treatment cycles. At the beginning of each cycle, we obtained the nanoparticle diameter that minimizes the tumor diameter before the next cycle starts. As the tumor size reduces, it became possible for multiple nanoparticle diameters to induce antitumor activity close to that of the minimizer of tumor diameter. For this reason, we opted to utilize the maximizer of tumor nanoparticle accumulation for the last treatment cycle, which maintained continuous tumor regression -similar to what could have happened by using the minimizer of tumor diameter -while gaining toxicological benefit. This protocol might be of use when treating tumors through multiple nanoparticle injections, i.e., minimizing the tumor diameter in the first few cycles until the size decreases to the point at which drug potency becomes less sensitive to nanoparticle diameter. Injection of nanoparticles with designs associated with high tumoral accumulation would then be optimal. Here, daily injection is assumed. The holiday time will be optimized in future studies.
On a patient-specific basis, the tumor model can be calibrated to fit morphological and phenotypical parameters of actual tumors after obtaining this information through imaging and histological analysis 17,[57][58][59][60][61][62] . Incorporation of further biological data, such as tumor vessel density, the amount of blood perfusion, integrity of endothelial cell layer, extracellular matrix protein, interstitial pressure, infiltrated immune cells, rate of tumor cell proliferation, resistant/anti-apoptotic mechanisms in tumor cells and their heterogeneities inside the tumor, would be expected to refine the prediction of drug delivery and therapeutic potency. If drug design is to be considered, genomic characterization of mutated cells could be incorporated in an expanded tumor model. Based on the predictive model, optimization studies, similar to those presented here, would then be conducted to find the optimal nanoparticle designs with clinically-based objectives, e.g., minimization of tumor size. Practically, the optimal nanoparticle design may not be clinically available, and thus the most similar currently available configuration would need to be chosen. In this case, this particular (sub-optimal) design would need to undergo modeling evaluation to obtain a revised prediction of potency. It may also prove useful to update the predictive model according to the actual tumor response, and perform repetitive optimization studies based on the updated model to design subsequent treatment cycles. This process could be repeated before each nanoparticle treatment, taking advantage of the speed of the optimization process, until tumor remission is achieved or neoadjuvant therapy is completed.
This study complements previous efforts aimed at characterizing nanoparticles that have maximal targeting properties. Sen Gupta 29 provided a review of computational and experimental investigations that were performed to find nanoparticle designs that enhance the margination, adhesion, and internalization process. While those studies explored designs based on evaluating finite possibilities of nanoparticle sizes and aspect ratios, finding an optimal design among a large number of candidates would only be possible through either expensive evaluation of all the possible designs or utilization of numerical optimization methods. The obtained optimal solutions in our study showed agreement with previous experimental findings. For instance, Joshi et al. 9 evaluated in vivo the tumoral uptake of a set of liposomes differing in size, mainly {60, 80, 200, 650, 670 nm}. Among the considered sizes, they found that the maximal uptake by tumors happened with 200 nm, a size that is closest to the maximizer (147-288 nm) of tumor nanoparticle accumulation obtained in this study. In 17 , van de Ven et al. studied the effect of drug loading and nanoparticle concentration on tumor regression, showing that there exists a nonlinear relationship between between tumor regression and drug concentration. In this study, we have incorporated the drug potency as an optimization variable to not only couple its selection with the nanoparticle design, but also to exploit its impact for minimizing the tradeoff between tumor reduction and nanoparticle accumulation. Charoenphol et al. 3 and Boso et al. 6 studied the effect of nanoparticle size on the binding to vasculature walls in vitro, postulating that maximum binding happens with a moderate value among the sizes tested. Our results shows that the relation between nanoparticle accumulation and size is non-monotonic, i.e., agreeing with 3,6 that the maximizer of nanoparticle binding is an interior optimum. From a therapeutic perspective, the minimizer of tumor diameter obtained in this study is 288 nm, which is very close to the 283 nm diameter nanoparticle synthesized with maximal entrapment efficiency in 8 . Not only is this size associated with a high drug load, but it is also suitable for passive tumor targeting through the EPR effect 8 . Rostami et al. 7 targeted tumors derived from the SKBR3 breast cancer cell line in mice. They used multiple injections of Doxorubicin-loaded nanoparticles, showing that the rate of tumor regression decreases with tumor size. Our study is consistent with this finding, showing that the change in tumor size decreases in time following multiple treatment cycles (Fig. 7). Having to dynamically select the proper nanoparticle size as the tumor changes in response to treatment requires different optimal nanoparticle designs at the beginning of each treatment cycle, which would be possible via optimization studies such as the ones presented here.
Cancer nanotherapy faces major barriers hindering its efficacy and limiting its translation to the clinic 63 . These barriers relate not only to synthesis processes and cost, but also to biological factors such as heterogeneity of the tumor microenvironment, immune system interactions, and acquired drug resistance 64 . Heterogeneity in therapeutic efficacy has been found among various tumor types, among patients bearing same tumor types, within single patients bearing multiple metastatic tumors or even within single tumors. There also exists time-and therapy-dependent evolution in tumor response to the same therapies. This study has applied optimization as a first step towards addressing these challenges to systematically help evaluate nanoparticle therapeutic potential by SCiEntiFiC REPORTs | (2018) 8:17768 | DOI:10.1038/s41598-018-35675-y maximizing drug targeting while reducing systemic toxicity. Ideally, these methods would allow for the design of nanoparticles to be customized to patient-specific tumors, while offering the capability to quickly evaluate adjustments to nanotherapy parameters during treatment. A major premise in this study is that nanoparticle-mediated drug delivery can offer substantially reduced systemic toxicity while increasing local drug delivery. Studies have shown that nanomedicine is associated with significant increases in the area under the drug concentration-time curve as compared to conventional chemotherapy 65 . This has spurred the development of nanomedications (around 9 FDA-approved and more than 30 under clinical trials 66 ) to treat various types of cancer. The analyses presented here are based on several assumptions. First, some biological parameters, such as hematocrit and drug diffusivity, were held constant. In reality, these parameters may vary by cancer type and tissue morphology and have an impact on the nanoparticle pharmacokinetics 37,67 . Future studies will use probability density functions to capture these variations and produce more robust designs that could be used for a wider range of tumors. Another source of uncertainty is variability in nanoparticle synthesis. Although optimal designs were generated in this study, it would be impractical to produce all nanoparticles equivalently, as they would be made with certain variations in their structure. Future work should account for these variations through robust optimization techniques that increase the applicability of the obtained design. In addition, tumor size at the beginning of the treatment may affect the optimal solution. Therefore, future work will study the impact of tumor size on the response. Finally, shape and functionalization have an effect on nanoparticle distribution and therapeutic potential 22,29 . This study shows that the integration of optimization in the design process makes it possible to investigate the efficiency of different nanoparticle designs, exploring new trends that may lower treatment toxicity while efficiently eradicating tumors.

Conclusion
In this study, the design of drug-carrying nanoparticles is optimized to maximize tumor regression and minimize the treatment toxicity. Tumor regression is quantified as the percentage change of tumor diameter to that at the beginning of the treatment. The treatment toxicity is measured as the fraction of the injected nanoparticles that accumulate in the tumor. The proposed nanoparticle designs provide a basis for further experimental and computational investigations. The study sheds light on design practices that increase nanoparticle circulation time while maintaining large drug encapsulation efficiency. This work lays a foundation to quantitatively evaluate preclinical nanoparticle-based drug delivery trials and support decisions in precision medicine where optimal solutions are required on a patient-specific basis. To achieve this goal, the following are suggested to focus the translational effort of mathematical modeling and optimization of nanotherapy. Models of adequate fidelity should be selected depending on tumor characteristics and the clinical problem under consideration. For instance, local injection of nanoparticles may benefit from designs that are suitable for treating tumors locally, for which the vasculature may be used for nanoparticle injections similar to intra-arterial chemotherapy for liver cancer 68 . Possible expansion of the model to accommodate a wider type of cancers include consideration of the extracellular matrix, presence of drug resistant cells, immune system activity, molecular-scale information (e.g., genetic, and proteomic), and nanoparticle accumulation in other organs such as the spleen, kidneys, and liver. A corresponding expansion of model input parameters would be necessary to capture this information. Lastly, a clinically-relevant user interface would enable parameter input and the delivery of modeling results.

Data Availability
All data analysed during this study are included in this published article. Additional datasets generated are available from the corresponding author upon request.