Pharmacokinetic—Pharmacodynamic Modeling of Tumor Targeted Drug Delivery Using Nano-Engineered Mesenchymal Stem Cells

Abstract

Nano-engineered mesenchymal stem cells (nano-MSCs) are promising targeted drug delivery platforms for treating solid tumors. MSCs engineered with paclitaxel (PTX) loaded poly(lactide-co-glycolide) (PLGA) nanoparticles (NPs) are efficacious in treating lung and ovarian tumors in mouse models. The quantitative description of pharmacokinetics (PK) and pharmacodynamics (PD) of nano-MSCs is crucial for optimizing their therapeutic efficacy and clinical translatability. However, successful translation of nano-MSCs is challenging due to their complex composition and physiological mechanisms regulating their pharmacokinetic-pharmacodynamic relationship (PK–PD). Therefore, in this study, a mechanism-based preclinical PK–PD model was developed to characterize the PK–PD relationship of nano-MSCs in orthotopic A549 human lung tumors in SCID Beige mice. The developed model leveraged literature information on diffusivity and permeability of PTX and PLGA NPs, PTX release from PLGA NPs, exocytosis of NPs from MSCs as well as PK and PD profiles of nano-MSCs from previous in vitro and in vivo studies. The developed PK–PD model closely captured the reported tumor growth in animals receiving no treatment, PTX solution, PTX-PLGA NPs and nano-MSCs. Model simulations suggest that increasing the dosage of nano-MSCs and/or reducing the rate of PTX-PLGA NPs exocytosis from MSCs could result in improved anti-tumor efficacy in preclinical settings.

1. Introduction

Since their discovery by A.J. Friedenstein in early 1970s [1], mesenchymal stem cells (MSCs) have been widely used in various therapeutic areas [2], especially in regenerative medicine [3,4]. Based on unique characteristics such as high growth and differentiation potential [5,6,7] and low immunogenicity [8,9], MSCs have been investigated in over 700 clinical trials worldwide for various clinical applications (https://clinicaltrials.gov/) [10]. MSCs exhibit specific tumor tropism in response to inflammatory signals secreted by tumor-infiltrated neutrophils and macrophages [11]. The C-X-C chemokine receptor type 4 (CXCR4)/stromal-derived factor-1 (SDF-1) axis is also involved in the recruitment of MSCs to the tumor site [12]. These mechanisms have been explored for tumor-targeted drug delivery using MSCs [13,14,15,16]. Early studies with MSCs typically involved genetic modifications to express anti–tumor proteins [17,18,19]. More recently, MSCs have been synthetically modified with polymeric nanoparticles (nano-MSCs) to enable the loading and delivery of small molecule anticancer drugs [15,16].

Unlike other targeted therapies, nano-MSCs take advantage of the inherent tumor homing capability of MSCs, to selectively and actively deliver potent chemotherapeutic drugs to tumor sites [2]. However, overexpression of drug efflux transporters such as P-glycoprotein prevents the direct and efficient loading of soluble forms of chemotherapeutic drugs, such as paclitaxel (PTX) [2,16,20]. In order to overcome this drawback, PTX is first encapsulated into biodegradable and biocompatible poly(lactide–co–glycolide) nanoparticles (PLGA NPs). MSCs are then engineered by simple incubation with PTX-loaded PLGA NPs (PTX-PLGA NPs) to create nano-MSCs (Figure 1). Nanoengineering does not affect the morphology, size and/or phenotype of MSCs [16]. Following intravenous administration, nano-MSCs home in to tumors within 15 min of the injection and their accumulation increases over time [15,16]. PTX is released slowly from nanoparticles and there is minimal drug release during the manufacturing process (~4 h), suggesting the stability of nanoparticles inside the cells during their preparation and tumor homing [15,16]. While a fraction of loaded PTX could release inside MSCs and result in cytotoxicity, our study showed that loading PTX in MSCs results in only 8–10% cell death [16]. MSCs are resistant to PTX-induced cytotoxicity, likely because MSCs overexpress P-glycoprotein, which actively effluxes PTX out of MSCs [16]. Our group has also shown that nano-MSCs exhibit greater anti-tumor efficacy and reduced toxicity at substantially lower PTX doses, compared to that with PTX solution and PTX-PLGA NPs (9 mg/kg total dose of PTX for nano-MSCs vs. 120 mg/kg total dose of PTX for free- and nanoparticle-encapsulated drugs) [15]. However, the pharmacokinetic (PK)/pharmacodynamic (PD) relationships and the mechanism behind the superior anti-tumor efficacy of nano-MSCs have not been fully characterized.

A quantitative understanding of the dose-exposure-response relationship, via mechanism-based models, is a crucial step to optimize the preclinical efficacy of developmental therapeutics. Developed with integrated data from in vitro and in vivo studies as well as literature-derived parameter values, mechanism-based PK–PD models are able to extrapolate drug exposure and therapeutic effects, which is anticipated to reduce experimental attrition and facilitate the translation of complex therapeutics [21,22,23]. Numerous mechanistic PK–PD models have been developed over the last decade for therapies using small molecule compounds, monoclonal antibodies and antibody drug conjugates [24,25,26]. However, very few studies have focused on modeling cell-based therapies [27,28].

Considering there are different components responsible for the delivery of PTX to tumor cells using nano-MSCs (Figure 1), it becomes crucial to understand and characterize the biodistribution of each of these components (PTX in the form of free PTX, PLGA NPs and nano-MSCs) to achieve an optimized dosing regimen with maximum therapeutic potential. However, it is difficult to fully understand the complex processes involved in the cellular and physiological disposition of MSCs, nanoparticles and PTX, without the use of a mechanism-based PK–PD model. Therefore, the aim of this study was to develop a mechanism-based PK–PD model to characterize the time course of PTX in the form of free drug, PLGA NPs and nano-MSCs and the tumor growth profiles following nano-MSC treatment in an A549 orthotopic mouse model of lung carcinoma. The developed model will be valuable for investigating the disposition of nano-MSCs and maximizing their therapeutic potential to facilitate the translation of nano-MSCs from bench to bedside.

2. Materials and Methods

2.1. Data

PK (tissue drug concentrations) and PD (tumor growth over time) data were obtained from our previously published study [15], which investigated the preclinical biodistribution and efficacy of nano-MSCs in an orthotopic model of A549 lung carcinoma. The in vivo studies in mice were performed according to the protocols approval by Institutional Animal Care and Use Committee (IACUC), University of Minnesota (Protocol ID: 1605-33821A). The lung tumors were developed in female Fox Chase SCID^® beige mice (Charles River Laboratories) through tail vein injection of A549-luc cells and tumor growth was monitored using in vivo bioluminescence imaging as described below.

2.1.1. PK Data

The tumor-bearing mice were randomized into 3 treatment groups, receiving a single 5 µg (0.25 mg/kg) intravenous (IV) injection of PTX solution, PTX-PLGA NPs, and nano-MSCs. Following the dose, mice (n = 3) were sacrificed at 2 h, and at 1-, 2-, 5- and 12-days post-treatment and blood, lung (tumor-bearing site), liver and spleen tissue samples were obtained. PTX concentrations were measured in these samples by LC/MS/MS as described before [15]. Concentrations measured as 0 in the original PK study were excluded during the model development process.

2.1.2. PD Data

Tumor growth was monitored by measuring tumor-associated bioluminescence using the IVIS Spectrum in vivo Imaging System (Caliper Life Sciences, Hopkinton, MA, USA). Treatments were initialized when tumor size reached 10⁵–10⁶ photon/sec. The tumor-bearing mice were randomized into 6 treatment groups receiving saline, untreated MSCs, MSCs with drug free blank NPs, PTX solution, PTX-PLGA NPs, or nano-MSCs. Mice in the PTX solution and PTX-PLGA NPs groups received 40 mg/kg of PTX on days 0, 4, and 8. Mice in the untreated MSCs, MSCs with blank NPs and nano-MSCs treatment groups received 1 × 10⁶ MSCs (for nano-MSCs: equivalent to 25 µg or 1.25 mg/kg PTX) on day 0 and 0.5 × 10⁶ MSCs every 14 days. For modeling purposes, animals receiving saline, untreated MSCs and MSCs loaded with blank NPs were pooled as the treatment-free group and used as control, given that PTX is the active tumor-killing cytotoxic drug.

2.2. PK–PD Modeling

Nonlinear mixed-effect modeling (NONMEM^® version 7.4 software; ICON Development Solutions, Ellicott City, MD, USA) was performed for all analyses using the first-order (FO) and Laplacian methods for the development of the PK and PK–PD models, respectively. Since the PK study was conducted with a destructive sampling approach (one sample per subject), a naïve pooled approach was utilized. Thus, no between-subject variability (BSV) terms were estimated in the PK models, and all the variability was attributed to the residual unexplained variability (RUV). The PTX concentrations measured as ng/g were converted to ng/mL by assuming a 1 g/mL tissue density [29]. A sequential PK–PD modeling approach was employed, in which PK models for the PTX solution, PTX-PLGA NPs and nano-MSCs were developed. The estimated PK parameters were then placed into NONMEM control stream and used to drive the tumor killing in the PD model. For the BSV in PK–PD model, a log-normal distribution was assumed. For the RUV in either PK or PK–PD models, proportional error models were used and assumed to be normally distributed with mean 0 and variance σ². Final NONMEM control streams for all analyses are provided in the Supporting Information. Exploratory analyses and diagnostic graphics were performed using R 3.6.3 (The R Foundation for Statistical Computing) and RStudio 1.1.453 (RStudio, Inc., Boston, MA, USA). Visual prediction checks (VPCs) and sampling importance resampling (SIR) were performed with Perl-speaks-NONMEM (PsN 4.9.0, Uppsala University, Uppsala, Sweden) under the Pirana^® interface [30].

2.2.1. PK Models

Following the administration of nano-MSCs, the PK of PTX is complicated by the simultaneous existence of PTX in forms of free PTX, PLGA NPs and nano-MSCs (Figure 1). In order to account for PTX co-existing in three different forms, a PK model was developed using 3 parallel layers, wherein each layer corresponds to each form of PTX. A “layer by layer” approach was applied as shown in Figure 2. A free PTX–PK model (bottom layer) was first developed using PK data of PTX solution treatment group. The estimated PK parameters were then fixed and a PTX-PLGA NPs layer (middle layer) was then grafted upon. PK data of PTX-PLGA NPs was then used to develop a PK model for PTX-PLGA NPs. The estimated PK parameters for the middle layer were then fixed and used to develop the nano-MSCs (top layer) PK model. The PK parameters in top layer were then estimated using the PK data of nano-MSCs.

PK Model for PTX Solution (Bottom Layer)

Figure 2A describes the PK model for the PTX solution. The PTX solution PK profiles in plasma and lung were represented by a three-compartment model. The model included a central compartment, a peripheral compartment and a tumor compartment with linear elimination from the central compartment (Eli). The model was used to simultaneously characterize the disposition of PTX in the systemic circulation and tumor-bearing lungs. The distribution of PTX between the central and peripheral compartments was represented by an inter-compartmental distribution process (Dist). The model characterized the disposition of PTX from the central compartment to the tumor compartment via vascular exchange (VE) and surface exchange (SE) processes as described previously by Shah et al. [31,32]. While the surface exchange describes the diffusion of therapeutic modalities across the tumor surface, the vascular exchange describes the permeation of therapeutic modalities across the blood vessel endothelium. The rates of these two processes were influenced by the size of therapeutic modalities and the size of tumor [33,34,35]. Due to the relatively similar scale of molecular weight of small molecule compounds, the rate constants of diffusion and permeability (D_PTX, P_PTX) and tumor accessible fraction (E_PTX) were assumed to be the same for all small molecules as shown in Table 1. The detailed PK model equations of PTX solution are provided below (Equations (1)–(3)).

$$\begin{gathered} \frac{dPT{X}_{central}}{dt}=CL{D}_{PTX}\times \left(\frac{A_{PT{X}_{peripheral}}}{V_{PT{X}_{peripheral}}}-\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right)-C{L}_{PTX}\times \left(\frac{A_{PT{X}_{central}} }{V_{PT{X}_{central}}}\right)-\frac{2\times P_{PTX}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \\ \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)-\frac{6\times D_{PTX}}{R_{tumor}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right) \end{gathered}$$

(1)

$$I{C}_{PT{X}_{central}}=Dos{e}_{PTX} IC:Initial Condition$$

$$\begin{array}{cc}\hfill \frac{dPT{X}_{tumor}}{dt}& =\frac{2\times P_{PTX}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)\hfill \\ & +\frac{6\times D_{PTX}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)\hfill \end{array}$$

(2)

$$I{C}_{PT{X}_{tumor}}=0$$

$$\frac{dPT{X}_{peripheral}}{dt}=-CL{D}_{PTX}\times \left(\frac{A_{PT{X}_{peripheral}}}{V_{PT{X}_{peripheral}}}-\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right)$$

(3)

$$I{C}_{PT{X}_{peripheral}}=0$$

where, PTX_central, PTX_tumor and PTX_peripheral represent central, tumor and peripheral compartment for PTX free drug, respectively. A_PTXcentral, A_PTXtumor and A_{PTXperipheral} describe the PTX amounts in central, tumor and peripheral compartment, respectively. The definitions of all the other parameters are shown in Table 1.

PK Model for PTX -PLGA NPs (Middle Layer)

Figure 2B describes the integrated PK model for PTX in the form of PLGA NPs. The model accounts for both PTX free drug and PTX in the form of PLGA NPs using the similar structural model that was previously used for PTX solution. The model consists of two layers, with each layer described using a three-compartment model with a central compartment, a peripheral compartment and a tumor compartment with linear elimination from the central compartment. The two layers of model are connected by a first order drug release process. The PTX release rate constant (K_rel) was derived from a separate in vitro study conducted previously [16], as described in the section K_rel and K_exo Estimation below. Detailed PK model equations for PTX in the form of PLGA NPs are provided below (Equations (4)–(9)).

$$\begin{array}{cc}\hfill \frac{dN{P}_{central}}{dt}& =CL{D}_{NP}\times \left(\frac{A_{N{P}_{peripheral}}}{V_{N{P}_{peripheral}}}-\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\right)-C{L}_{NP}\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\right) \hfill \\ & -\frac{2\times P_{NP}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)\hfill \\ & -\frac{6\times D_{NP}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)-K_{rel}\times A_{N{P}_{central}}\hfill \end{array}$$

(4)

$$I{C}_{N{P}_{central}}=Dos{e}_{PTXNP}$$

$$\begin{array}{cc}\hfill \frac{dN{P}_{tumor}}{dt}& =\frac{2\times P_{NP}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)\hfill \\ & +\frac{6\times D_{NP}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)-K_{rel}\times A_{N{P}_{tumor}} \hfill \end{array}$$

(5)

$$I{C}_{N{P}_{tumor}}=0$$

$$\frac{dN{P}_{peripheral}}{dt}=-CL{D}_{NP}\times \left(\frac{A_{N{P}_{peripheral}}}{V_{N{P}_{peripheral}}}-\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\right)-K_{rel}\times A_{N{P}_{peripheral}}$$

(6)

$$I{C}_{N{P}_{peripheral}}=0$$

$$\begin{array}{cc}\hfill \frac{dPT{X}_{central}}{dt}& =CL{D}_{PTX}\times \left(\frac{A_{PT{X}_{peripheral}}}{V_{PT{X}_{peripheral}}}-\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right)-C{L}_{PTX}\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right) \hfill \\ & -\frac{2\times P_{PTX}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)\hfill \\ & -\frac{6\times D_{PTX}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)+K_{rel}\times A_{N{P}_{central}}\hfill \end{array}$$

(7)

$$I{C}_{PT{X}_{central}}=0$$

$$\begin{array}{cc}\hfill \frac{dPT{X}_{tumor}}{dt}& =\frac{2\times P_{PTX}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)\hfill \\ & +\frac{6\times D_{PTX}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)+K_{rel}\times A_{N{P}_{tumor}} \hfill \end{array}$$

(8)

$$I{C}_{PT{X}_{tumor}}=0$$

$$\frac{dPT{X}_{peripheral}}{dt}=-CL{D}_{PTX}\times \left(\frac{A_{PT{X}_{peripheral}}}{V_{PT{X}_{peripheral}}}-\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right)+K_{rel}\times A_{N{P}_{peripheral}}$$

(9)

$$I{C}_{PT{X}_{peripheral}}=0$$

where, NP_central, NP_tumor, NP_peripheral, PTX_central, PTX_tumor and PTX_peripheral represent the central, tumor and peripheral compartments for PTX in the form of PLGA NPs and central, tumor and peripheral compartments for PTX free drug, respectively in PTX-PLGA NPs PK model. A_NPcentral, A_NPtumor, A_NPperipheral, A_PTXcentral, A_PTXtumor and A_{PTXperipheral} describe the PTX amounts in central, tumor and peripheral compartments for PTX in the form of PLGA NPs and for PTX free drug, respectively. The definitions of all the other parameters are shown in Table 1.

PK Model for Nano-MSCs (Top Layer)

Figure 2C describes the integrated PK model for PTX in the form of nano-MSCs. The model accounts for PTX free drug, PTX in the form of PLGA NPs and PTX in the form of nano-MSCs using the similar structural model previously used for PTX solution and PTX-PLGA NPs. The model consists of three layers. The structures of the PTX free drug and PTX-PLGA NPs layers are the same as the PK model for PTX-PLGA NPs. The nano-MSCs layer is described using three compartments: a central compartment, a peripheral compartment and a tumor compartment. Two unidirectional extravasation processes connect the central compartment with the peripheral and the tumor compartment in nano-MSCs layer, assuming the extravasation of nano-MSCs follows single direction. The nano-MSC and PTX-PLGA-NP layers of model are connected via a first order exocytotic process. The PLGA NPs’ exocytotic rate constant (K_exo) was derived from a separate in vitro study conducted previously [36] as described in the section K_rel and K_exo Estimation below. The rate constant of PTX free drug release from nano-MSCs is assumed to be the same as the PTX free drug release from PLGA NPs layer and described via a first order drug release rate constant, K_rel. The assumption is based on the relatively fast efflux of PTX free drug from MSCs. Detailed PK model equations for PTX in the form of nano-MSCs are provided below (Equations (10)–(18)).

$$\frac{dMS{C}_{central}}{dt}=-K_{ct}\times A_{MS{C}_{central}}-K_{cp}\times A_{MS{C}_{central}}-K_{rel}\times A_{MS{C}_{central}}-K_{exo}\times A_{MS{C}_{central}}$$

(10)

$$I{C}_{MS{C}_{central}}=Dos{e}_{nanoMSCs}$$

$$\frac{dMS{C}_{tumor}}{dt}=K_{ct}\times A_{MS{C}_{central}}-K_{rel}\times A_{MS{C}_{tumor}}-K_{exo}\times A_{MS{C}_{tumor}}$$

(11)

$$I{C}_{MS{C}_{tumor}}=0$$

$$\frac{dMS{C}_{peripheral}}{dt}=K_{cp}\times A_{MS{C}_{central}} -K_{rel}\times A_{MS{C}_{peripheral}}-K_{exo}\times A_{MS{C}_{peripheral}}$$

(12)

$$I{C}_{MS{C}_{peripheral}}=0$$

$$\begin{array}{cc}\hfill \frac{dN{P}_{central}}{dt}& =CL{D}_{NP}\times \left(\frac{A_{N{P}_{peripheral}}}{V_{N{P}_{peripheral}}}-\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\right)-C{L}_{NP}\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\right)\hfill \\ & -\frac{2\times P_{NP}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)\hfill \\ & -\frac{6\times D_{NP}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)-K_{rel}\times A_{N{P}_{central}}+K_{exo}\times A_{MS{C}_{central}} \hfill \end{array}$$

(13)

$$I{C}_{N{P}_{central}}=0$$

$$\begin{array}{cc}\hfill \frac{dN{P}_{tumor}}{dt}& =\frac{2\times P_{NP}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)\hfill \\ & +\frac{6\times D_{NP}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{N{P}_{central} }}{V_{N{P}_{central}}}\times E_{NP}-\frac{A_{N{P}_{tumor}}}{VT}\right)-K_{rel}\times A_{N{P}_{tumor}}+K_{exo}\times A_{MS{C}_{tumor}} \hfill \end{array}$$

(14)

$$I{C}_{NPtumor}=0$$

$$\frac{dN{P}_{peripheral}}{dt}=-CL{D}_{NP}\times \left(\frac{A_{N{P}_{peripheral}}}{V_{N{P}_{peripheral}}}-\frac{A_{N{P}_{central}}}{V_{N{P}_{central}}}\right)-K_{rel}\times A_{N{P}_{peripheral}}+K_{exo}\times A_{MS{C}_{peripheral}}$$

(15)

$$I{C}_{N{P}_{peripheral}}=0$$

$$\begin{array}{cc}\hfill \frac{dPT{X}_{central}}{dt}& =CL{D}_{PTX}\times \left(\frac{A_{PT{X}_{peripheral}}}{V_{PT{X}_{peripheral}}}-\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right)-C{L}_{PTX}\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right)\hfill \\ & -\frac{2\times P_{PTX}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)\hfill \\ & -\frac{6\times D_{PTX}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)+K_{rel}\times A_{N{P}_{central}}+K_{rel}\times A_{MS{C}_{central}}\hfill \end{array}$$

(16)

$$I{C}_{PT{X}_{central}}=0$$

$$\begin{array}{cc}\hfill \frac{dPT{X}_{tumor}}{dt}& =\frac{2\times P_{PTX}\times R_{Cap}}{R_{Krogh}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)\hfill \\ & +\frac{6\times D_{PTX}}{R_{Tumor}{}^2}\times VT\times \left(\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\times E_{PTX}-\frac{A_{PT{X}_{tumor}}}{VT}\right)+K_{rel}\times A_{N{P}_{tumor}}+K_{rel}\times A_{MS{C}_{tumor}} \hfill \end{array}$$

(17)

$$I{C}_{PT{X}_{tumor}}=0$$

$$\frac{dPT{X}_{peripheral}}{dt}=-CL{D}_{PTX}\times \left(\frac{A_{PT{X}_{peripheral}}}{V_{PT{X}_{peripheral}}}-\frac{A_{PT{X}_{central}}}{V_{PT{X}_{central}}}\right)+K_{rel}\times A_{N{P}_{peripheral}}+K_{rel}\times A_{MS{C}_{peripheral}}$$

(18)

$$I{C}_{PT{X}_{peripheral}}=0$$

In the above equations, MSC_central, MSC_tumor and MSC_peripheral represent the central, tumor and peripheral compartments, respectively, for PTX in the form of nano-MSCs. NP_central, NP_tumor, NP_peripheral, PTX_central, PTX_tumor and PTX_peripheral represent the central, tumor and peripheral compartments for PTX in the form of PLGA NPs and for PTX free drug, respectively. A_MSCcentral, A_MSCtumor and A_{MSCperipheral}, A_Npcentral, A_Nptumor and A_Npperipheral, A_PTXcentral, A_PTXtumor and A_{PTX peripheral} describe the PTX amounts in central, tumor and peripheral compartments for PTX in the form of nano-MSCs, PLGA NPs and PTX free drug, respectively. The definitions of all the other model parameters are presented in Table 1.

K_rel and K_exo Estimation

To estimate K_rel, free PTX release from PLGA NPs was assumed to follow a first order release process and the in vivo PTX release kinetics from PLGA NPs was assumed to be comparable to in vitro PTX release kinetics. The in vitro PTX release from PLGA NPs study conducted by Sadhukha et al. [16] was used to derive K_rel. A first order association model was fitted using least square method with GraphPad Prism 7. Additionally, the efflux of free PTX from nano-MSCs was assumed to happen instantaneously and therefore drug release from PLGA NPs within the cells was considered the rate-limiting step for free PTX release from nano-MSCs. Thus, the derived K_rel was used for describing free drug release from both PLGA NPs and nano-MSCs.

To estimate K_exo, exocytosis of PLGA NPs from nano-MSCs was assumed to follow a first order exocytosis process. Moreover, the in vivo PLGA NP exocytosis was assumed to be comparable to in vitro exocytosis. The in vitro NPs exocytosis study conducted by Layek et al. [36] was used to derive the rate constant for NP exocytosis (K_exo). A first order decay model was fitted using least square method with GraphPad Prism 7.

2.2.2. PK–PD Model

The estimated PK parameters from free PTX, PTX-PLGA NPs and nano-MSCs–PKmodels were fixed and made available in the NONMEM control stream for the tumor growth model. The predicted PTX concentrations for free drug, PLGA NPs and nano-MSCs were assumed to inhibit tumor growth as shown in Figure 2D. Tumor inhibition caused by PTX solution and PTX-PLGA NPs was described using Emax model (Equations (19) and (20)). IC50 values for PTX solution and PTX-PLGA NPs were previously reported by in vitro cytotoxicity studies [16]. Tumor inhibition caused by PTX in the form of nano-MSCs (Equation (21)) was described using a linear model. The tumor growth was represented as an exponential tumor growth function along with the tumor inhibition caused by PTX free drug, PTX-PLGA NPs and nano-MSCs. Detailed tumor growth PD model equations are provided below (Equations (19)–(22)).

$$Kkil{l}_{PTX}= \frac{Kma{x}_{PTX}\times \left(\frac{A_{PT{X}_{tumor}}}{VT}\right)}{IC{50}_{PTX}+\frac{A_{PT{X}_{tumor}}}{VT}}$$

(19)

$$Kkil{l}_{NP}= \frac{Kma{x}_{NP}\times \left(\frac{A_{N{P}_{tumor}}}{VT}\right)}{IC{50}_{NP}+\frac{A_{N{P}_{tumor}}}{VT}}$$

(20)

$$Kkil{l}_{MSC}= Kma{x}_{MSC}\times \left(\frac{A_{MS{C}_{tumor}}}{VT}\right)$$

(21)

$$\frac{dTV}{dt}=Kg0\times TV-Kkil{l}_{PTX}\times TV-Kkil{l}_{NP}\times TV-Kkil{l}_{MSC}\times TV$$

(22)

$$I{C}_{TV}=TVBL$$

In the above equations, Kkill_PTX, Kkill_NP and Kkill_MSC are the tumor killing rates caused by PTX free drug, PTX in the form of PLGA NPs and PTX in the form of nano-MSCs, respectively. VT represents the tumor volume and TV represents the tumor compartment with the unit consistent with tumor bioluminescence. The conversion between TV and VT was derived from a separate animal study as described in the section Correlation between Tumor Weight and Tumor Bioluminescence. The definitions of all the other PD model parameters are provided in Table 2.

The tumor bioluminescence profiles from animals receiving no treatment, PTX solution, PTX–PLGA NPs and nano-MSCs were fitted simultaneously with a population model. In order to account for the potentially early dropouts of animals having a high tumor burden, a modified M3 method-based approach proposed by Martin et al. was utilized [37]. In this approach, the highest tumor bioluminescence measurement from each animal was set as the upper “detection” limit for that animal, while measurements below the upper limit were treated as continuous data. Tumor bioluminescence measurements that were not available due to early dropouts, were assumed to be collected as regular sampling scheme and assumed to be above the upper limit for each animal. These measurements were treated as categorical data to maximize the likelihood with respect to model parameters. This method was proposed to overcome the underestimation of treatment effects introduced by the different dropout time of animals in different treatment arms.

Correlation between Tumor Weight and Tumor Bioluminescence

The animal studies were conducted according to the protocols approval by Institutional Animal Care and Use Committee (IACUC) at the University of Minnesota (Protocol ID: 1605-33821A). Six to eight-week-old female Fox Chase SCID Beige mice were purchased from Charles River Laboratories. The orthotopic human lung cancer model was developed in these mice by the injection of A549-luc cells (1 × 10⁶ cells) through the intravenous (tail vein) route. The tumor bioluminescence (photon/sec) and tumor weights (g) were reported for A549 tumor mice (n = 8) before and after euthanasia, respectively. Assuming the density of tumor tissue as 1 g/mL, the relation between tumor bioluminescence in photon/sec and tumor volume in mL can be investigated. In order to prevent the existence of negative tumor volume during the PD modeling process, a power function (Equation (23)) was used to investigate the association between tumor volume and tumor bioluminescence with Microsoft Excel.

$$VT=a\times T{V}^{b }$$

(23)

where, VT is tumor volume in mL; and TV is tumor bioluminescence in 10⁶ photon/s.

2.3. Model Evaluations

PK and PK–PD models were evaluated based on the objective function values (OFVs), the stability of the model, and goodness-of-fit plots. The precisions of the final model parameters were evaluated using relative standard error (RSE) generated from covariance steps and sampling-importance-resampling (SIR)-based 95% confidence intervals (CIs) [38,39]. The final PK–PD models were evaluated using prediction-corrected visual predictive checks (pcVPC; 1000 simulations) [40].

2.4. PK and PK–PD Model Simulations

After the development of PK model for nano-MSCs, simulations were performed following the administration of 5 µg (0.25 mg/kg) PTX equivalent nano-MSCs as a single IV bolus dose into the central compartment of the nano-MSCs layer. PK profiles of the total PTX, PTX in the form of free drug, PLGA NPs and nano-MSCs were simulated in both plasma and lung. In addition, following the development of the full PK–PD model, simulations (n = 100) were carried out for animals receiving nano-MSCs, using typical values and BSV estimated in the final models. Different dosing scenarios were simulated to reflect: (i) different dosing intervals of the same total dose (10⁶ MSCs (equivalent to 25 µg or 1.25 mg/kg PTX) on day 0 and 0.5 × 10⁶ MSCs every two weeks; 10⁶ MSCs on day 0 and 0.25 × 10⁶ MSCs every week; 10⁶ MSCs on day 0 and 0.125 × 10⁶ MSCs every 3.5 days); (ii) different doses with the same dosing interval (10⁶ MSCs on day 0 and 0.25 × 10⁶ MSCs every two weeks; 10⁶ MSCs on day 0 and 0.5 × 10⁶ MSCs every two weeks; 10⁶ MSCs on day 0 and 0.75 × 10⁶ MSCs every two weeks; 10⁶ MSCs on day 0 and 2 × 10⁶ MSCs every two weeks). The effect of different K_rel (0.00425, 0.0085, 0.017 h⁻¹) and K_exo (0.06, 0.081, 0.1 h⁻¹) values on tumor growth were also simulated with the original dosing regimen of 10⁶ MSCs on day 0 and 0.5 × 10⁶ MSCs every two weeks. All the simulations were performed with mrgsolve package in R. The median and 10th and 90th percentiles of simulated tumor bioluminescence were calculated and plotted using R.

3. Results

3.1. PK and PD Data Exploration

The drug exposure in the lung (tumor-bearing tissue) following a single dose of nano-MSCs is much higher compared to that with PTX solution and PTX-PLGA NPs with comparable plasma drug exposure (Figure 3A). The tumor growth profiles are shown in Figure 3B. The relatively uneven width of 95% confidence intervals at each time point indicates relatively noisy measurements. It is noted that tumor growth profiles of animals receiving no treatments exhibit exponential growth up to around 63 days (1512 h). A plateau was observed for animals receiving no treatments after 63 days because of the death of animals with relatively high tumor bioluminescence. Similar plateaus were observed for animals in the other treatment groups.

3.2. In Vitro and In Vivo Kinetics Parameter Estimation

Given the limited PK data, the rate constant K_rel and K_exo were estimated with data acquired from in vitro drug release and NPs exocytosis studies [16,36]. The first order association model used estimated K_rel as 0.0085 h⁻¹ with R² of 0.845 (Figure 4A), whereas the first order decay model used estimated K_exo as 0.56 h⁻¹ with R² of 0.734 (Figure 4B). Since PK models were developed with respect to PTX mass, the estimated K_exo was further multiplied by the reported PTX loading of PTX-PLGA NPs, 0.148 mg PTX/mg NP [16]. The resulted rate constant 0.083 h⁻¹ was fixed in the following nano-MSCs–PK model development process. However, a fixed K_exo of 0.083 h⁻¹ was found to induce relatively imprecise estimates (large RSE%) of the other model parameters. Hence, we tested different K_exo values ranging from 0.080 to 0.085 h⁻¹ and found that fixing K_exo to 0.081 h⁻¹ provided a reasonable RSE% for all model parameters and gave the lowest objective function value (OFV). The association between tumor weight and tumor bioluminescence is shown in Figure 4C.

3.3. PK and PK–PD Model Parameters

The total PTX concentrations in plasma and lung following the administration of PTX solution, PTX-PLGA NPs and nano-MSCs were reasonably predicted using the developed PK models in both plasma and tumor-bearing lungs (Figure 5). The PK model parameter estimations, their RSE, SIR-based 95% CIs and sources of literature-derived parameters for PTX solution, PTX-PLGA NPs and nano-MSCs are provided in Table 1. The final PK–PD model parameter estimates, their RSE, SIR-based 95% CIs, and sources of literature-derived parameters for animals receiving no treatment, PTX solution, PTX-PLGA NPs and nano-MSCs are provided in Table 2.

3.4. Model Evaluations

The overall goodness-of-fit plots indicated little evidence to reject the PK–PD model (Supplementary Figures S1–S5). The prediction-corrected visual prediction check (pcVPC) was stratified by treatment groups (no treatment, treatments of PTX solution, PTX-PLGA NPs and nano-MSCs) as shown in Figure 6 on a log scale. The pcVPC suggests the developed PK–PD model is able to predict the observed tumor bioluminescence reasonably well. The final PK models and PK–PD model parameter estimates were all within their SIR-based 95% CIs (Table 2).

Table 1. Parameters for pharmacokinetic model.

Parameters	Estimates (%RSE)	SIR Medians (95% CIs)	Units	Sources	Definitions
Model Associated Parameters
PPTX	0.0875		cm/h	Literature [31,33]	Permeability rate constant for PTX free drug
DPTX	0.01		cm2/h	Literature [31,33]	Diffusion rate constant for PTX free drug
EPTX	0.44		unitless	Literature [31,33]	Tumor fraction accessible by PTX free drug
Rkrogh	0.0008		cm	Assumed [41]	Inter-capillary distance
Rcap	0.0075		cm	Assumed [42]	Radius of tumor-associated capillaries
VT	0.3 (PK); Dynamic (PD)		mL	Assumed (PK); Dynamic (PD)	Tumor volume
Rtumor	0.42 (PK); Dynamic (PD)		cm	Calculated/Assumed (PK); Dynamic (PD)	Tumor radius
CLPTX	0.909 (23%)	0.952 (0.684, 1.305)	mL/h	Estimated	Clearance for PTX free drug
CLDPTX	0.336 (68%)	0.424 (0.193, 0.712)	mL/h	Estimated	Distribution clearance for PTX free drug
VPTXcentral	6.64 (30%)	7.06 (2.81, 11.49)	mL	Estimated	Central compartment volume of distribution for PTX free drug
VPTXperipheral	18.5 (53%)	22.6 (12.6, 35.7)	mL	Estimated	Peripheral compartment volume of distribution for PTX free drug
fuPTX	0.0237 (46%)	0.0237 (0.0119, 0.0499)	unitless	Estimated	Plasma to blood ratio for PTX free drug
PNP	0.00035		cm/h	Literature [43]	Permeability rate constant for PTX in the form of PLGA NPs
DNP	3.6 × 10−6		cm2/h	Literature [44]	Diffusion rate constant for PTX in the form of PLGA NPs
ENP	0.055		unitless	Literature [33]	Tumor fraction accessible by PTX in the form of PLGA NPs
Krel	0.0085		1/h	Calculated	First order drug release rate constant for PTX free drug from PTX-PLGA NPs
CLNP	0.241 (17%)	0.259 (0.172, 0.340)	mL/h	Estimated	Clearance for PTX in the form of PLGA NPs
CLDNP	0.0627 (48%)	0.0681 (0.0349, 0.1244)	mL/h	Estimated	Distribution clearance for PTX in the form of PLGA NPs
VNPcentral	1.32 (17%)	1.41 (0.85, 1.96)	mL	Estimated	Central compartment volume of distribution for PTX in the form of PLGA NPs
VNPperipheral	43.2 (141%)	79.0 (16.6, 219.3)	mL	Estimated	Peripheral compartment volume of distribution for PTX in the form of PLGA NPs
fuNP	0.00302 (57%)	0.00325 (0.00127, 0.00684)	unitless	Estimated	Plasma to blood ratio for PTX in the form of PLGA NPs
Kexo	0.081		1/h	Calculated	First order exocytosis rate constant for PTX-PLGA NPs from nano-MSCs
Kct	1.45 (22%)	1.39 (1.04, 1.84)	1/h	Estimated	Rate constant describing central to tumor compartment transfer for PTX in the form of nano-MSCs
Kcp	10.2 (2%)	10.2 (9.8, 10.8)	1/h	Estimated	Rate constant describing central to peripheral compartment transfer for PTX in the form of nano-MSCs
VMSCcentral	7.15 × 10−8 (30%)	8.50 × 10−8 (2.22 × 10−8, 1.94 × 10−7)	mL	Estimated	Central compartment volume of distribution for PTX in the form of nano-MSCs
VMSCperipheral	15021 (55%)	13565 (932, 33,564)	mL	Estimated	Peripheral compartment volume of distribution for PTX in the form of nano-MSCs
Residual Unexplained Variability (RUV, proportional, CV%)
εPTX_plasma	116% (20%)	118.5% (67.7%, 176.6%)		Estimated	RUV for PTX solution plasma PK profiles
εPTX_tumor	64.9% (25%)	71.5% (46.4%, 103.9%)		Estimated	RUV for PTX solution lung PK profiles
εNP_plasma	81.2% (71%)	100.8% (50.0%, 167.6%)		Estimated	RUV for PTX-PLGA NPs plasma PK profiles
εNP_tumor	54.5% (14%)	58.8% (42.7%, 74.0%)		Estimated	RUV for PTX-PLGA NPs lung PK profiles
εMSC_plasma	87.3% (26%)	86.4% (73.4%, 118.5%)		Estimated	RUV for nano-MSCs plasma PK profiles
εMSC_tumor	58.1% (22%)	61.2% (43.4%, 85.7%)		Estimated	RUV for nano-MSCs lung PK profiles

Table 1. Parameters for pharmacokinetic model.

Parameters	Estimates (%RSE)	SIR Medians (95% CIs)	Units	Sources	Definitions
Model Associated Parameters
PPTX	0.0875		cm/h	Literature [31,33]	Permeability rate constant for PTX free drug
DPTX	0.01		cm2/h	Literature [31,33]	Diffusion rate constant for PTX free drug
EPTX	0.44		unitless	Literature [31,33]	Tumor fraction accessible by PTX free drug
Rkrogh	0.0008		cm	Assumed [41]	Inter-capillary distance
Rcap	0.0075		cm	Assumed [42]	Radius of tumor-associated capillaries
VT	0.3 (PK); Dynamic (PD)		mL	Assumed (PK); Dynamic (PD)	Tumor volume
Rtumor	0.42 (PK); Dynamic (PD)		cm	Calculated/Assumed (PK); Dynamic (PD)	Tumor radius
CLPTX	0.909 (23%)	0.952 (0.684, 1.305)	mL/h	Estimated	Clearance for PTX free drug
CLDPTX	0.336 (68%)	0.424 (0.193, 0.712)	mL/h	Estimated	Distribution clearance for PTX free drug
VPTXcentral	6.64 (30%)	7.06 (2.81, 11.49)	mL	Estimated	Central compartment volume of distribution for PTX free drug
VPTXperipheral	18.5 (53%)	22.6 (12.6, 35.7)	mL	Estimated	Peripheral compartment volume of distribution for PTX free drug
fuPTX	0.0237 (46%)	0.0237 (0.0119, 0.0499)	unitless	Estimated	Plasma to blood ratio for PTX free drug
PNP	0.00035		cm/h	Literature [43]	Permeability rate constant for PTX in the form of PLGA NPs
DNP	3.6 × 10−6		cm2/h	Literature [44]	Diffusion rate constant for PTX in the form of PLGA NPs
ENP	0.055		unitless	Literature [33]	Tumor fraction accessible by PTX in the form of PLGA NPs
Krel	0.0085		1/h	Calculated	First order drug release rate constant for PTX free drug from PTX-PLGA NPs
CLNP	0.241 (17%)	0.259 (0.172, 0.340)	mL/h	Estimated	Clearance for PTX in the form of PLGA NPs
CLDNP	0.0627 (48%)	0.0681 (0.0349, 0.1244)	mL/h	Estimated	Distribution clearance for PTX in the form of PLGA NPs
VNPcentral	1.32 (17%)	1.41 (0.85, 1.96)	mL	Estimated	Central compartment volume of distribution for PTX in the form of PLGA NPs
VNPperipheral	43.2 (141%)	79.0 (16.6, 219.3)	mL	Estimated	Peripheral compartment volume of distribution for PTX in the form of PLGA NPs
fuNP	0.00302 (57%)	0.00325 (0.00127, 0.00684)	unitless	Estimated	Plasma to blood ratio for PTX in the form of PLGA NPs
Kexo	0.081		1/h	Calculated	First order exocytosis rate constant for PTX-PLGA NPs from nano-MSCs
Kct	1.45 (22%)	1.39 (1.04, 1.84)	1/h	Estimated	Rate constant describing central to tumor compartment transfer for PTX in the form of nano-MSCs
Kcp	10.2 (2%)	10.2 (9.8, 10.8)	1/h	Estimated	Rate constant describing central to peripheral compartment transfer for PTX in the form of nano-MSCs
VMSCcentral	7.15 × 10−8 (30%)	8.50 × 10−8 (2.22 × 10−8, 1.94 × 10−7)	mL	Estimated	Central compartment volume of distribution for PTX in the form of nano-MSCs
VMSCperipheral	15021 (55%)	13565 (932, 33,564)	mL	Estimated	Peripheral compartment volume of distribution for PTX in the form of nano-MSCs
Residual Unexplained Variability (RUV, proportional, CV%)
εPTX_plasma	116% (20%)	118.5% (67.7%, 176.6%)		Estimated	RUV for PTX solution plasma PK profiles
εPTX_tumor	64.9% (25%)	71.5% (46.4%, 103.9%)		Estimated	RUV for PTX solution lung PK profiles
εNP_plasma	81.2% (71%)	100.8% (50.0%, 167.6%)		Estimated	RUV for PTX-PLGA NPs plasma PK profiles
εNP_tumor	54.5% (14%)	58.8% (42.7%, 74.0%)		Estimated	RUV for PTX-PLGA NPs lung PK profiles
εMSC_plasma	87.3% (26%)	86.4% (73.4%, 118.5%)		Estimated	RUV for nano-MSCs plasma PK profiles
εMSC_tumor	58.1% (22%)	61.2% (43.4%, 85.7%)		Estimated	RUV for nano-MSCs lung PK profiles

Table 2. Parameters for pharmacodynamic model.

Parameters	Estimates (%RSE)	SIR Medians (95% CIs)	Units	Sources	Definitions
Model Associated Parameters
Kg0CTR	0.00339 (6%)	0.00339 (0.00325, 0.00352)	/h	Estimated	First order tumor growth rate constant for animals receiving no treatments
Kg0PTX	0.00372 (9%)	0.00367 (0.00338, 0.00397)	/h	Estimated	First order tumor growth rate constant for animals receiving PTX solution
Kg0PTXNP	0.00417(5%)	0.00418 (0.00394, 0.00445)	/h	Estimated	First order tumor growth rate constant for animals receiving PTX-PLGA NPs
Kg0MSC	0.00588 (11%)	0.00585 (0.00528, 0.00639)	/h	Estimated	First order tumor growth rate constant for animals receiving nano-MSCs
TVBLCTR	0.360 (16%)	0.365 (0.251, 0.488)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving no treatments
TVBLPTX	0.376 (43%)	0.387 (0.208, 0.593)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving PTX solution
TVBLPTXNP	0.539 (18%)	0.549 (0.365, 0.730)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving PTX-PLGA NPs
TVBLMSC	0.227 (14%)	0.232 (0.173, 0.287)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving nano-MSCs
KmaxPTX	0.00343 (26%)	0.00332 (0.00261, 0.00398)	/h	Estimated	Maximal tumor killing rate induced by PTX free drug
KmaxNP	0.000427 (372%)	0.000755 (0.000035, 0.002192)	/h	Estimated	Maximal tumor killing rate induced by PTX in the form of PLGA NPs
KMSC	4.35 × 10−6 (30%)	4.31 × 10−6 (2.79 × 10−6, 5.75 × 10−6)	/(h*(ng/mL))	Estimated	Linear tumor killing rate constant induced by PTX in the form of nano-MSCs
IC50PTX	1.5		ng/mL	Literature [16]	Concentration of PTX free drug can introduce 50% KmaxPTX
IC50NP	5.7		ng/mL	Literature [16]	Concentration of PTX in the form of PLGA NPs can introduce 50% KmaxNP
Between Subject Variability (BSV, proportional, CV%)
ηTVBL	96.4% (9%)	99.1% (77.3%, 123.6%)	Shrinkage (0%)	Estimated	BSV on baseline tumor bioluminescence
Residual Unexplained Variability (RUV, CV%)
εPTX_plasma	49.3% (5%)	49.5% (45.1%, 53.7%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving no treatments
εPTX_tumor	66.8% (7%)	67.4% (59.1%, 76.6%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving PTX solution
εNP_plasma	62.4% (8%)	63.0% (54.9%, 71.2%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving PTX–PLGA- NPs
εNP_tumor	62.2% (10%)	63.0% (54.8%, 72.2%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving nano-MSCs

Note: * indicate multiplication.

Table 2. Parameters for pharmacodynamic model.

Parameters	Estimates (%RSE)	SIR Medians (95% CIs)	Units	Sources	Definitions
Model Associated Parameters
Kg0CTR	0.00339 (6%)	0.00339 (0.00325, 0.00352)	/h	Estimated	First order tumor growth rate constant for animals receiving no treatments
Kg0PTX	0.00372 (9%)	0.00367 (0.00338, 0.00397)	/h	Estimated	First order tumor growth rate constant for animals receiving PTX solution
Kg0PTXNP	0.00417(5%)	0.00418 (0.00394, 0.00445)	/h	Estimated	First order tumor growth rate constant for animals receiving PTX-PLGA NPs
Kg0MSC	0.00588 (11%)	0.00585 (0.00528, 0.00639)	/h	Estimated	First order tumor growth rate constant for animals receiving nano-MSCs
TVBLCTR	0.360 (16%)	0.365 (0.251, 0.488)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving no treatments
TVBLPTX	0.376 (43%)	0.387 (0.208, 0.593)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving PTX solution
TVBLPTXNP	0.539 (18%)	0.549 (0.365, 0.730)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving PTX-PLGA NPs
TVBLMSC	0.227 (14%)	0.232 (0.173, 0.287)	106 photon/s	Estimated	Baseline tumor bioluminescence for animals receiving nano-MSCs
KmaxPTX	0.00343 (26%)	0.00332 (0.00261, 0.00398)	/h	Estimated	Maximal tumor killing rate induced by PTX free drug
KmaxNP	0.000427 (372%)	0.000755 (0.000035, 0.002192)	/h	Estimated	Maximal tumor killing rate induced by PTX in the form of PLGA NPs
KMSC	4.35 × 10−6 (30%)	4.31 × 10−6 (2.79 × 10−6, 5.75 × 10−6)	/(h*(ng/mL))	Estimated	Linear tumor killing rate constant induced by PTX in the form of nano-MSCs
IC50PTX	1.5		ng/mL	Literature [16]	Concentration of PTX free drug can introduce 50% KmaxPTX
IC50NP	5.7		ng/mL	Literature [16]	Concentration of PTX in the form of PLGA NPs can introduce 50% KmaxNP
Between Subject Variability (BSV, proportional, CV%)
ηTVBL	96.4% (9%)	99.1% (77.3%, 123.6%)	Shrinkage (0%)	Estimated	BSV on baseline tumor bioluminescence
Residual Unexplained Variability (RUV, CV%)
εPTX_plasma	49.3% (5%)	49.5% (45.1%, 53.7%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving no treatments
εPTX_tumor	66.8% (7%)	67.4% (59.1%, 76.6%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving PTX solution
εNP_plasma	62.4% (8%)	63.0% (54.9%, 71.2%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving PTX–PLGA- NPs
εNP_tumor	62.2% (10%)	63.0% (54.8%, 72.2%)		Estimated	RUV for tumor bioluminescence profiles in animals receiving nano-MSCs

Note: * indicate multiplication.

3.5. Nano-MSC PK Model Simulation

Simulations were performed to examine the anticipated PK profiles following the administration of nano-MSCs. Following the administration of 5 µg PTX equivalent nano-MSCs, total PTX concentration, PTX concentrations in the form of PTX solution, PLGA NPs and nano-MSCs in both plasma and lung were simulated (Figure 7). Our simulations show that following a single dose of nano-MSCs, PTX in the form of nano-MSCs quickly depleted from plasma and rapidly accumulated in the tumor-bearing lung. In the tumor-bearing lung, PTX in the form of nano-MSCs showed a mono-exponential decline. In addition, 64 h following the nano-MSC administration, free PTX starts to dominate the PTX concentration in the tumor-bearing lung.

3.6. PK–PD Model Simulation

The effect of nano-MSCs administration on tumor growth during a relevant dose regimen was examined using a simulation approach based on our PK–PD model. Our simulation suggested that different dosing intervals have little impact on the overall effectiveness of nano-MSCs (Figure 8A). On the other hand, simulation with different doses showed a significant difference in tumor growth inhibition, with the greatest inhibition driven by the highest dosage (Figure 8B). We further evaluated the impact of K_rel (Figure 8C) and K_exo (Figure 8D) of nano-MSCs on tumor growth. Our simulation suggested that K_rel has a little impact while K_exo has a significant impact on tumor growth profile. A smaller K_exo is associated with a greater tumor inhibition.

4. Discussion

The present study characterized the time course of free PTX, PTX-PLGA NPs and nano-MSCs and tumor growth profiles following nano-MSC dosing in an orthotopic mouse model of lung carcinoma using a nonlinear mixed-effects modeling approach. A mechanism-based PK–PD model for nano-MSCs was developed by leveraging data collected from previous in vitro and in vivo studies [15,16,36]. The developed preclinical PK models were able to adequately predict PTX exposure in animals following the administration of PTX solution, PTX-PLGA NPs and nano-MSCs in both plasma and the tumor-bearing lung. The developed preclinical PK models for PTX solution, PTX-PLGA NPs and nano-MSCs was merged with a tumor growth model to develop a PK–PD model that can characterize tumor bioluminescence–time profiles in animals receiving no treatment, PTX solution, PTX-PLGA NPs and nano-MSCs.

Mechanistic models that account for the disposition of targeted anticancer agents in the tumor-bearing tissue are often more suitable for developing a translatable PK–PD relationship. Xie et al. [45] have reported several difficulties in building a PK–PD model with plasma drug exposure only, which points to the necessity of collecting drug concentration in tumor-bearing tissues and characterizing the relationship between drug exposure at the tumor site and tumor growth. In this study, we have developed one such model that described the disposition of PTX in the form of free drug, PTX-PLGA NPs and nano-MSCs in both plasma and tumor-bearing lungs following nano-MSC dosing. We used a destructive sampling approach to measure PTX concentrations at various tissue sites. In this approach, one animal was sacrificed per time point as previously described [15]. We used a naïve pooled modeling approach to develop PK models for PTX solution, PTX-PLGA NPs and nano-MSCs [46]. This study design allowed us to characterize the drug exposure in both plasma and tumor sites, which facilitated the development of a mechanistic PK–PD model.

Following the administration of nano-MSCs, the PK of PTX is complicated by the simultaneous existence of PTX in forms of free PTX, PLGA NPs and nano-MSCs (Figure 1). In order to account for PTX co-existing in three different forms, a PK model was developed by merging 3 parallel PK models (layers), wherein each layer correspond with each form of PTX. The total PTX concentrations in plasma and lung following the administration of PTX solution, PTX-PLGA NPs and nano-MSCs were reasonably predicted using the developed PK models in both plasma and tumor-bearing lungs (Figure 5). In order to implement the “layer by layer” approach, several assumptions were made to describe PTX release and PLGA-NP exocytosis from nano-MSCs. PTX release from PLGA NPs and nano-MSCs were assumed to follow a first order release process. This is considered a limitation of the study since it may over-simplify the physiological conditions. A more physiologically relevant drug release kinetic model such as Higuchi model may better characterize the drug release profiles from matrix-based system such as PLGA NPs and decrease parameter imprecision [47,48].

Once the PK model was developed, a mechanistic PK–PD model was further developed to characterize the tumor growth in animals receiving no treatment, PTX solution, PTX-PLGA NPs and nano-MSCs, using a population modeling approach. Previous studies have frequently modeled tumor growth empirically using a mix of first order and zero order growth functions [32,49,50]. When we tried to implement similar models, convergence issues were repeatedly reported and model parameters associated with the linear growth function were estimated with insufficient precision. When we further explored tumor bioluminescence data in control (no treatment) group, we found that only few data points reflect linear growth (Figure 3B, red trace; Supplemental Figure S1, animal ID 1-18). The scarcity of data acquired in the linear (zero-order) growth phase may reflect the difference in the tumor growth pattern in different preclinical tumor models. Our model was developed using data collected from orthotopic lung tumor model rather than subcutaneous models. Orthotopic models are considered more clinically relevant than subcutaneous models [51,52]. However, the potential high aggressiveness of orthotopic tumor models [53] may cause most of animals to die before tumor volume reaches the linear growth phase. As a result, an exponential (first-order) growth function was only used to characterize tumor growth in our model.

In preclinical studies, subcutaneous tumor models are often used to study tumor growth for the purpose of optimizing the therapeutic benefits of newly developed cancer therapies. In these studies, calipers are frequently used to monitor tumor growths [31,54,55], whereas in orthotopic models, non-invasive bioluminescence measurements are needed to monitor tumor growth. Bioluminescence measurements are generally considered as a semi-quantitative method for measuring tumor size [45]. The accuracy of measurements can be compromised by numerous intrinsic and extrinsic factors such as the administration routes of substrate and cellular environments [45,56]. As a result, bioluminescence measurements are typically noisy (Figure 3B) and difficult to describe using a tumor growth PD model. This was noticeable during the development of our PK–PD model. We experienced difficulties in achieving model convergence. This is partially explained by excessive experimental noise introduced by bioluminescence measurements. The inclusion of a variability term on parameters such as exponential growth rate constant (Kg0) led to unsuccessful convergence of the model. As shown in Table 2, only one between subject variability (BSV) term on the baseline tumor volume was included in the final model. The failure to estimate multiple BSV terms in the final model may also reflect the difficulties in modeling tumor bioluminescence measurements.

One limitation of the PK–PD model was that a few parameters related to PTX-PLGA NPs such as V_NPperipheral and Kmax_NP were estimated with relatively poor precision as shown in Table 1 and Table 2. The unexpected large relative standard error (RSE) may potentially reflect an unoptimized sampling scheme, insufficient data for model development and/or over-parameterization of the structural model. This is due to the fact that the original study was designed for the purpose of demonstrating the superior targeting and efficacy achieved by nano-MSCs and not for the purpose of PK–PD modeling. Further studies are needed to provide additional experimental PK–PD data to better characterize PTX–PLGA-NP model parameters with high precision.

Given the complex composition of nano-MSCs, examining the PK profiles of different PTX components of this therapy experimentally is extremely difficult. Our simulations performed using the developed nano-MSCs–PK model provided mechanistic insights about the disposition of the three different forms of PTX as well as total PTX, following nano-MSCs administration. In tumor-bearing mice, following nano-MSCs administration, the relatively fast targeting, retention and extravasation of MSCs quickly deplete nano-MSCs from plasma to peripheral tissues, especially tumor-bearing sites. Nano-MSCs deliver large amount of PTX to tumor-bearing organ despite premature release of free PTX and PTX-PLGA NPs in systemic circulation (Figure 7). Subsequently, free drug and PTX-PLGA NPs released from nano-MSCs at the tumor-bearing site gradually diffuse back into the systemic circulation, which explains the shallow ascending phase following the initial rapid decline reflected in total PTX plasma concentrations after nano-MSCs dosing (Figure 3A and Figure 7A). While it would be valuable to determine the concentrations of PTX in various forms (free-PTX, PTX-PLGA NPs and PTX-MSCs), it is extremely difficult to accurately estimate the various fractions because of practical and analytical limitations. For example, separating free drug from nanoparticles requires either dialysis or ultracentrifugation. Dialysis is time consuming (nanoparticles will continue to release the drug during dialysis) while ultracentrifugation-based methods are confounded by presence of blood cells and other proteins. Because of these difficulties, most studies involving nano drug carriers report total drug concentration [57]. As more advanced analytical and separation technologies become available, future studies could analyze the concentrations of each of the forms separately to gain a better understanding of the disposition of nano-MSCs and to validate our model predictions.

Simulations performed using the developed PK–PD model can be useful in optimizing preclinical dosing regimens for nano-MSCs. Previous studies have suggested that cancer therapies, such as gefitinib and trametinib, can achieve a better therapeutic efficacy in treating certain types of cancer when administered in multiple low doses compared with a single high dose [58,59]. However, our simulations suggest that reducing the dosing interval of nano-MSC treatment while maintaining the same total dose does not improve efficacy. In contrast, increasing the total dose of nano-MSCs while maintaining the same dosing intervals is anticipated to improve the therapeutic efficacy of nano-MSCs (Figure 8A,B). In the case of nano-MSCs, only 9 mg/kg of PTX total dose was delivered over the course of the study, which is much lower than clinically administered doses for PTX. Because our current dose is low, we do expect to be able to inject higher amounts of PTX with nano-MSCs with minimal toxicity. However, the maintenance dose of nano-MSCs is limited by the relatively large size of this therapeutic modality [60]. From our experimental experience, administrating an IV dose greater than 4 × 10⁶ nano-MSCs may increase the chance of forming cell clots within blood vessels and cause immediate death of experimental animals [15,61]. Thus, a maximal maintenance dose of 2 × 10⁶ MSCs (equivalent to 50 μg or 2.5 mg/kg PTX) was simulated. It has to be noted that if animals can tolerate such multiple doses, it would result in complete tumor remission. However, the goal should be to improve the payload capacity of nano-MSCs so as to minimize dosing frequency to further improve their effectiveness. We believe that efficacy studies conducted with the simulated maintenance doses are needed to evaluate the practicability of the predicted dosing regimens and to validate the model predictions. We expanded our simulations of the PK–PD model to examine properties that may optimize the therapeutic potential of nano-MSCs. Achieving more sustained release of chemotherapeutic drugs could improve their therapeutic efficacy [62,63]. However, in our studies, the rate of drug release from PLGA NPs had little impact on the overall tumor growth inhibition (Figure 8C). On the other hand, our model suggests that reducing the rate of exocytosis of PLGA NPs from nano-MSCs may improve nano-MSC efficacy (Figure 8D).

5. Conclusions

The mechanistic model developed closely captures the complex PK–PD behavior of nano-MSCs in an orthotopic A549 lung carcinoma model. The developed model provides an opportunity to examine the effect of various parameters to further optimize the therapeutic potential of nano-MSCs. The model also provides mechanistic insights regarding the in vivo kinetics of nano-MSCs and can help guide future preclinical study design. The modeling workflow illustrated by the development of the current PK–PD model can serve as a platform for developing a mechanistic model for nano-MSCs in other tumor models. This will further help to optimize the therapeutic efficacy and clinical translatability of nano-MSCs.