A Hybrid In Silico and Tumor-on-a-Chip Approach to Model Targeted Protein Behavior in 3D Microenvironments

Simple Summary Engineered proteins possess a great therapeutic potential, but the development of such therapies is impeded during preclinical studies by the lack of in vitro models that accurately simulate the human physiology. Animal models, on the other hand, also have difficulties predicting human responses, and are ethically concerning. In this study, we employed a hybrid approach where we combined mathematical modeling with 3D in vitro models that mimic aspects of the tumor microenvironment, in order to simulate the delivery of therapeutic proteins targeting cancer cells and to predict the biological activity. By cross-comparing simulated and experimental data from 3D models, we were able to correctly predict the best dose needed to deliver toxic proteins specifically to tumor cells, while leaving the surrounding non-tumor cells untouched. This study shows the potential of combining computational approaches with novel in vitro models to advance the development of protein therapeutics. Abstract To rationally improve targeted drug delivery to tumor cells, new methods combining in silico and physiologically relevant in vitro models are needed. This study combines mathematical modeling with 3D in vitro co-culture models to study the delivery of engineered proteins, called designed ankyrin repeat proteins (DARPins), in biomimetic tumor microenvironments containing fibroblasts and tumor cells overexpressing epithelial cell adhesion molecule (EpCAM) or human epithelial growth factor receptor (HER2). In multicellular tumor spheroids, we observed strong binding-site barriers in combination with low apparent diffusion coefficients of 1 µm2·s−1 and 2 µm2 ·s−1 for EpCAM- and HER2-binding DARPin, respectively. Contrasting this, in a tumor-on-a-chip model for investigating delivery in real-time, transport was characterized by hindered diffusion as a consequence of the lower local tumor cell density. Finally, simulations of the diffusion of an EpCAM-targeting DARPin fused to a fragment of Pseudomonas aeruginosa exotoxin A, which specifically kills tumor cells while leaving fibroblasts untouched, correctly predicted the need for concentrations of 10 nM or higher for extensive tumor cell killing on-chip, whereas in 2D models picomolar concentrations were sufficient. These results illustrate the power of combining in vitro models with mathematical modeling to study and predict the protein activity in complex 3D models.


Introduction
Engineered proteins including antibodies and non-immunoglobulin binding proteins that target cancer-specific features play an increasingly important role in the development of novel cancer therapies [1]. Compared with small molecules, proteins offer enhanced specificity and/or selective recognition of surface features [2]. Furthermore, their small size offers benefits with respect to tissue accessibility compared with most nanoparticulate delivery systems [3]. Nevertheless, the efficient delivery of therapeutic proteins into tumors can be hindered because of a binding-site barrier [4], a high interstitial pressure as a result of poor tumor vasculature, a defective lymphatic system, and a dense extracellular matrix [5]. When the local concentration of the therapeutic agent is sublethal deep inside the tumor tissue, not only is a less efficacious treatment observed, but the development of therapy resistance is also much more likely [6,7].
Traditionally, tumor targeting has been studied using 2D culture models and animal models. However, the aforementioned barriers to efficient drug delivery are not replicated in 2D culture systems, and in vivo models suffer from a low throughput. Furthermore, in vivo models are difficult to control, expensive, less suitable for real-time uptake studies, and ethically problematic. In addition, there are increasing worries regarding the poor extrapolation of results from animal model outcomes to the human setting [8][9][10]. Consequently, there has been growing interest in the utilization of in vitro complex 3D tissues that mimic relevant aspects of human tumor microenvironments for studying drug delivery [11,12].
Multicellular tumor spheroids (MCTS) consisting of a cluster of tumor cells are a suitable 3D model to study drug delivery to avascular tumor regions [12]. MCTS are one of the simplest 3D models that can replicate tumor features, such as gradients of oxygen and waste products and the high density of the extracellular matrix [12,13]. MCTS can be generated with a high throughput and have been used to efficiently study the penetration of nanoparticles, antibodies, and antibody fragments in 3D [14][15][16]. Alternative models that include more complex features of the tumor microenvironment have gained more attention in recent years. In particular, microfluidic cancer-on-a-chip platforms, either consisting solely of cancer cells or co-cultured with stromal cells [17][18][19][20][21], provide a physiologically more relevant 3D microenvironment because of the ability to compartmentalize and better control gradients [22,23]. Furthermore, cancer-on-a-chip systems also allow for real-time imaging and quantitative analyses of the biological processes in ways unachievable with current animal models. Therefore, these models have the potential for a more in-depth investigation of the transport processes of drugs, including therapeutic proteins, into and around tumor tissues, and allow for a quick evaluation of novel formats with respect to their transport behavior. Importantly, the highly controlled nature of experiments with microfluidic models facilitates the application of modeling approaches for the description of transport processes, which, so far, have been limited to the description of processes in either static spheroids [14] or in vivo models [24][25][26][27][28].
In this study, we investigated the specific delivery of designed ankyrin repeat proteins (DARPins) directed against human epidermal growth factor receptor 2 (HER2) or epithelial cell adhesion molecule (EpCAM) in complex 3D tumor microenvironments. We included optically cleared spheroids and a microfluidic tumor model for real-time visualization of the tumor penetration, using 2D models as a reference. Based on our experimental setup, we also developed a mathematical model for evaluating protein delivery in spheroids and on-chip ( Figure 1). We used DARPins because they are highly stable engineered proteins that can be selected to bind, with very high affinities, to virtually any target [29]. Additionally, DARPins are 10× smaller than conventional antibodies, giving them better tissue-penetration properties [30]. To investigate the specific activity towards tumor cells in our tumor-on-a-chip system, we fused DARPins to the Pseudomonas aeruginosa exotoxin A (ETA) and compared the experimental data with mathematical predictions of an adapted mathematical model considering the cell-killing activity of the agents. The predictions of the model aligned with our experimental observations regarding the expected amount of cell death throughout the microfluidic chip after 48 h, illustrating the potential of a hybrid approach that combines tumor-on-a-chip systems with mathematical models for investigating the transport processes and biological activities of candidate protein therapies. The insert shows a close-up view of the microchannels that connect the side channels of the tumor compartment. (B) Spheroids were modelled in COMSOL as spheres into which the free DARPins (Da) diffused from outside to inside, while encountering receptors on the surface of the cells (Re). The model was simplified into a 1D model using the spherical expression of Equation (1) (see Material and Methods). B refers to DAR-Pin-receptor complex. (C) The tumor-on-a-chip 3D geometry was reproduced in COMSOL. DARPins diffuse from the side channel to the tumor compartment, which is connected by small microchannels. The matrix properties and cell density were taken into account by implementing distinct diffusion coefficients (abbreviated as Dwater and Dmatrigel) and available volume (void) fractions (ε) for the tumor compartment and for the side channel.

Cell Lines and Culture Media
The experiments were performed with human lung adenocarcinoma Calu-3 cells, human mammary gland ductal carcinoma BT-474 cells (both ATCC, Manassas, VA, USA), and two primary human fibroblast cell lines-IMR-90 (ATCC, Manassas, VA, USA ) and C5120-which were kindly donated by Dr. Rodenburg (Radboud University Medical Center, Nijmegen, The Netherlands) [31]. The BT-474 cells were cultured in Gibco Roswell Park Memorial Institute medium (RPMI), Calu-3 cells in a minimum essential medium Eagle (MEM), and the fibroblasts in Medium 199 (Gibco, Invitrogen, Carlsbad, CA, USA). The insert shows a close-up view of the microchannels that connect the side channels of the tumor compartment. (B) Spheroids were modelled in COMSOL as spheres into which the free DARPins (Da) diffused from outside to inside, while encountering receptors on the surface of the cells (Re). The model was simplified into a 1D model using the spherical expression of Equation (1) (see Material and Methods). B refers to DARPinreceptor complex. (C) The tumor-on-a-chip 3D geometry was reproduced in COMSOL. DARPins diffuse from the side channel to the tumor compartment, which is connected by small microchannels. The matrix properties and cell density were taken into account by implementing distinct diffusion coefficients (abbreviated as D water and D matrigel ) and available volume (void) fractions (ε) for the tumor compartment and for the side channel.

Cell Lines and Culture Media
The experiments were performed with human lung adenocarcinoma Calu-3 cells, human mammary gland ductal carcinoma BT-474 cells (both ATCC, Manassas, VA, USA), and two primary human fibroblast cell lines-IMR-90 (ATCC, Manassas, VA, USA ) and C5120-which were kindly donated by Dr. Rodenburg (Radboud University Medical Center, Nijmegen, The Netherlands) [31]. The BT-474 cells were cultured in Gibco Roswell Park Memorial Institute medium (RPMI), Calu-3 cells in a minimum essential medium Eagle (MEM), and the fibroblasts in Medium 199 (Gibco, Invitrogen, Carlsbad, CA, USA). All of the media were supplemented with 10% fetal calf serum (FCS). For the experiments with microfluidic devices, 1% v/v penicillin-streptomycin (Sigma-Aldrich, St. Louis, MO, USA)

Cloning, Expression and Conjugation of DARPins and DARPin Fusions
The EpCAM-binding DARPin Ec1 and the control DARPin Off7 with an N-terminal MRGS-His6 tag and a C-terminal cysteine and FLAG tag were expressed from pQIq vectors in the Escherichia coli (E. coli) strain BLR(DE3), as described previously [32]. The HER2-binding DARPins 9_26 was similarly cloned in pQIq through BamHI and HindIII restriction sites and was expressed in BLR(DE3). The DARPins fused to the translocation and catalytic domains of Pseudomonas aeruginosa exotoxin A (ETA; aa 252-612) were generated by replacing the consensus DARPin NI 3 C in Ec1-ETA (252-412)-NI 3 C-KDEL with the catalytic domain of ETA through SpeI and PstI restriction sites [33]. The EpCAM-binding DARPin Ec1 was subsequently replaced by the HER2-targeting DARPins 9_26, or the control DARPin Off7 through the BamHI and HindIII restriction sites. The DARPin-fusion proteins were expressed in the E. coli strain BL21(DE3) essentially, as described previously [33], with the exception that a 2YT medium (16 g·L −1 peptone, 10 g·L −1 yeast extract, 5 g·L −1 NaCl) containing 50 µg mL −1 ampicillin (Sigma-Aldrich, St. Louis, MO, USA) was used for the protein expression. DARPins were coupled to Alexa Fluor 680 (Thermo Fisher Scientific, Waltham, MA, USA), as described previously [32]. The purified DARPins and DARPin-toxin fusions are shown in Figure S1A,B, and a schematic of the engineered DARPins is presented in Figure 2A, as well as the confirmation of the specific binding to Calu-3 cells in Figure 2B.

Quantification of Receptor Density
The cells were detached using a 10 mM EDTA solution (Merck, Darmstadt, Germany) at 4 • C and were stained with a Phycoerythrin (PE)-labeled anti-EpCAM antibody (Abcam, Cambridge, MA, USA) or PE-labeled anti-HER2 antibody (Abcam, Cambridge, MA, USA). A PE-labeled IgG1 isotype control antibody (Life Technologies Carlsbad, CA, USA) was used as a negative control. The fluorescence intensity for each antibody was measured on a MACSQuant Analyzer 10 flow cytometer (Miltenyi, Bergisch Gladbach, Germany), and a calibration curve was generated using the BD QuantiBrite™ PE fluorescence quantitation kit (BD Biosciences, Franklin Lakes, NJ, USA), according to the manufacturer's instructions. The data were analyzed with FlowJo 10.0.8rl (BD Life Sciences, Ashland, OR, USA). For quantifying the receptor number ( Figure S2), we utilized the median of the whole cell population for EpCAM, whereas for the quantification of HER2 receptor density, we used the median of the positive cell population because of inhomogeneous staining. Through staining with receptor-targeting DARPins, we could confirm the homogeneous expression of HER2 on the cell lines (Calu-3: Figure 2C, BT-474: Figure S3), which led us to conclude that staining with the commercial antibodies was incomplete. Staining with DARPins also supported the similar number of binding sites found for EpCAM and HER2 on BT-474 cells ( Figure S3).

Determination of Binding of DARPin to Cell Surface Receptors
Calu-3, BT-474, and IMR-90 cells were seeded in 24-well culture plates. The cells were grown overnight and subsequently incubated for 2 h with Alexa Fluor 680-labeled anti-EpCAM or anti-HER2 DARPins at 37 • C. The cells were then washed twice with PBS and detached with PBS containing 5 mM EDTA for Calu-3 and BT-474. The IMR-90 cells required detachment with trypsin. The cells were centrifuged at 300× g for 3 min and resuspended in PBS. The cells were measured using the MACSQuant Analyzer 10 flow cytometer (Miltenyi, Bergisch Gladbach, Germany), and the data were analyzed with FlowJo 10.0.8rl (BD Life Sciences, Ashland, OR, USA). [33], with the exception that a 2YT medium (16 g·L peptone, 10 g·L yeast extract, 5 g·L NaCl) containing 50 μg mL −1 ampicillin (Sigma-Aldrich, St. Louis, MO, United States) was used for the protein expression. DARPins were coupled to Alexa Fluor 680 (Thermo Fisher Scientific, Waltham, MA, United States), as described previously [32]. The purified DARPins and DARPin-toxin fusions are shown in Figure S1A,B, and a schematic of the engineered DARPins is presented in Figure 2A, as well as the confirmation of the specific binding to Calu-3 cells in Figure 2B.

Specific DARPin Binding to Tumor Cells in Co-Cultures
Calu-3 cells and IMR-90 cells were stained with carboxyfluorescein succinimidyl ester (CFSE) or CellTrace yellow (both from Thermo Fisher Scientific, Waltham, MA, USA), respectively. The cells were seeded at a density of 50,000 cells/well for Calu-3 and 30,000 cells/well for IMR-90 in µ-slide eight-well chambers (Ibidi, Gräfelfing, Germany). Co-cultures were grown overnight and subsequently incubated for 2 h with various DARPin-Alexa Fluor 680 conjugates at 500 nM for 2 h. The cells were washed with complete medium (Medium 199, Life Technologies Carlsbad, CA, USA) and images were taken with an HC PL APO 63x/1.20 N.A. water immersion objective on a Leica TCS SP8 confocal microscope (Leica Microsystems, Wetzlar, Germany).

Toxicity of DARPin-Toxin Fusions
Calu-3 cells were seeded in a 96-well plate at a density of 20,000 cells/well. The cells were grown overnight and incubated for 48 h with DARPins fused to the translocation and catalytic domains of ETA. The cells were then either incubated with a 0.1 mg·mL −1 resazurin working solution in complete medium and incubated for 2 h, or grown for 96 h more before being incubated with resazurin. Afterwards, fluorescence was measured with 525-540 nm excitation and 620-640 nm emission light on a BioTek Synergy 2 Multi-Mode Reader (BioTek Instruments, Winooski, VT, USA). The relative viability was calculated from the fluorescent signals by subtracting the resazurin background signal and normalizing to untreated cells.

DARPin Penetration into Tumor Spheroids
BT-474 spheroids were produced via the hanging-drop method by pipetting drops of 30 µL with 5000 cells, as described previously [20]. Spheroids were allowed to form overnight at 37 • C in an incubator with 5% CO 2 . The DARPins were added in 1.5 µL to the 30 µL drops with the spheroid and were incubated for the indicated times at 37 • C. Spheroids were fixated for 2 h in 4% (w/v) PFA (Sigma-Aldrich, St. Louis, MO, USA) and kept in 2% PFA overnight at 4 • C. The spheroids were washed and embedded in a 20 µL drop of collagen solution (4 mg·mL −1 ) in an eight-well IBIDI chamber (Gräfelfing, Germany). The IBIDI chamber was kept at room temperature for 5 min, after which the collagen was polymerized for 40 min at 37 • C. A modified SeeDB protocol was employed for clearing, as described previously [20]. Confocal images were taken with an HCX PL APO 40x/0.85 dry objective or an HCX PL APO 10x/0.40 dry objective on a Leica TCS SP8 confocal microscope (Leica Microsystems, Wetzlar, Germany).

Quantification of DARPin Penetration Depth in Tumor Spheroids
A python-based macro for Fiji [34] was created and used to quantify the penetration depth of the fluorescently labeled DARPins into the BT-474 spheroids (Supplementary Methods). The slice representing the widest cross-section, and hence the maximum distance from spheroid rim to its center, was used for the quantification of penetration. First, the region of interest (ROI) was manually selected around the spheroid. After determining the ROI center, 360 radial lines covering 360 • were drawn from the center to the border of the ROI. Next, a Gaussian blur was applied before determining the signal intensity along the lines. The maximum intensity values determined for 360 lines were averaged. The half-maximal intensity was calculated, which was set as the rim of the spheroid, and was used to determine the distance of penetration. To determine the spheroid size, the radii were averaged. Tumor penetration was expressed as a percentage of the relative penetration from the outer rim of the spheroid.

Fabrication of Tumor-on-a-Chip
A design similar to the one reported by Chernyavska et al. was employed, composed of three parallel channels on glass on a device fabricated utilizing replica-molded polydimethylsiloxane (PDMS; Dow Corning, Midland, MI, USA) [20]. In this case, the spheroid traps were removed from the central chamber and 5 mm holes were punched to make the reservoirs connected to the side channels, using a biopsy puncher (Kai Industries, Seki, Gifu, Japan). All of the manufacturing steps were performed as reported previously [20].

Real-Time DARPin Penetration in the Tumor-on-a-Chip System
Co-cultures of IMR-90 or C5120 fibroblasts were labeled according to the manufacturer's instructions with CellTrace yellow, and CellTrace violet was used for the BT-474 cells. The labeled cells were then embedded in a 10 mg·mL −1 Matrigel (Corning, Corning, NY, USA) matrix and loaded manually in the tumor compartment of the microfluidic chip at a density of 2.0 × 10 7 cells mL −1 . The chips were placed in the incubator at 37 • C) for 15 min to allow the Matrigel to polymerize. After this time, a medium was added in the side channels and each of the four reservoirs were filled with the same volume of medium. We compared the Matrigel with collagen (4 mg·mL −1 ), but under the evaluated circumstances, we noted significant contraction of the collagen gel. After 48 h, protein diffusion through the Matrigel matrix and the specific targeting of cancer cells was investigated by adding Alexa Fluor 680-labeled DARPins in one of the side channels and filling up the two reservoirs with the exact same volume of medium immediately before imaging. Using the LASX software (Leica Microsystems, Wetzlar, Germany) with a Leica TCS SP8 confocal microscope (Leica Microsystems, Wetzlar, Germany) equipped with an HCX PL APO 10x/0.40 dry objective (Leica Microsystems, Wetzlar, Germany), one position per chip was selected and time lapses were recorded by taking an image of each marked position every 1 min for 1 h. For the image analyses and subsequent comparisons with the simulated data, Fiji was used. The Matrigel region of the middle channel was manually selected and plot profiles of fluorescence intensity across this channel were generated for the images at 0, 15, 30, 45, and 60 min. The intensity values were converted to concentrations using an averaged background-corrected maximum value in the side channel and the concentration at which DARPins were added (500 nM). We assumed a linear relationship between the fluorescence intensity and concentration.
For the toxicity studies in the tumor-on-a-chip, BT-474 cells labeled with CFSE and C5120 fibroblasts labeled with CellTrace violet were grown in the chip for 24 h. After this, DARPin-ETA fusions were added to one of the side channels of the chip at concentrations of 1 or 10 nM. The cell viability after 48 h of culture was determined by adding the cell death staining propidium iodide (PI; Sigma-Aldrich, St. Louis, MO, USA) into the chips and through quantification, as described below.

Quantification of Cell Death in Tumor-on-a-Chip System
To quantify the relative amount of cell death of each cell type in the tumor-on-a-chip after treatment with the DARPin-toxin fusions, a macro in Fiji was developed (Supplementary Methods). Briefly, each channel was filtered using a Gaussian blur. A threshold was applied to the images and the positive area corresponding to the cells or the PI staining was converted into a region of interest (ROI). The ROIs corresponding to one cell channel and PI were overlapped, and a ROI of the overlapping areas was created. Subsequently, the percentage of the overlapping area compared to the total cell area in the threshold image was calculated. The relative amount of cell death (based on surface area) was then reported as the ratio of the overlapping area between the PI channel and the cell channel.

Mathematical Model of DARPins Diffusion In Vitro
A binding-diffusion mathematical model was implemented by introducing the custom governing equations (Equations (1) to (3)) in the finite element software COMSOL Multiphysics v5.5 (Stockholm, Sweden) using the equation-based interface. A schematic of the input variables related to the therapy and the tissue, as well as the modelled outcomes, is given in Figure S4. The underlying coupled equations were initially reported by Thurber et al. [14] and Torres et al. [28], and describe the kinetics of unbound DARPin (1), receptor (2), and DARPin-receptor complex (3), respectively.
Where Da, Re, and B are the DARPin, receptor, and DARPin-receptor complex concentrations, respectively. The starting value of Da in the side channel reflects the chosen concentration and is indicated per experiment. The values for Re were experimentally determined or taken from the literature ( Table 1). The starting value of B is, by definition, zero.  D is the DARPin diffusion coefficient in ECM and was determined experimentally ( Figure S5A); ε is the interstitial volume fraction (see below for how values were acquired); and k on , k off , and k e are the binding, release, and endocytosis rate constants of the receptors, respectively and were taken from the literature (Table 1). Re 0 is the receptor concentration at t = 0. Note that the binding terms k on and k off must be divided by the available interstitial volume fraction for the interstitial diffusion, referred to as the void fraction (ε), because the DARPin-receptor interactions only occur in the space accessible to DARPins. ε is calculated as follows, ε = V_ECM/V_total, where V refers to volume and ECM to the extracellular matrix.
We considered the cells to be immobile, so that among the three variables, only Da can diffuse (D not equal to zero). The following two distinct assumptions could be made for the receptor regeneration term of R s : (i) If cells are not killed during DARPin diffusion, all receptors are recycled, R s = k e (Re + B). This assumption was taken when investigating the delivery of non-toxic DARPins. (ii) If cells are killed, internalized DARPin-receptor complexes are no longer recycled, R s = k e Re. This assumption was taken when investigating the delivery of fusion toxins, i.e., DARPins fused to a fragment of the bacterial toxin ETA.
Transient simulations were performed to predict the concentration profiles over time in both spheroids and tumor-on-chip systems.
Diffusion inside the spheroids was modelled on a spherical geometry, as shown in Figure 1A. The model was simplified into a 1D problem by taking the spherical coordinates' expression for the Laplace operator. For the boundary conditions, a constant surrounding the DARPin concentration was assumed at the sphere border, and zero-flux was assumed at the center (r = 0), so as to consider the spherical symmetry. The void fraction for modeling the delivery in spheroids was taken from the literature ( Table 1).
The tumor-on-chip model was solved on a 3D geometry, as shown in Figure 1B. To consider the presence of cells in the matrix, distinct values for the void fraction and diffusion coefficient (D water and D ECM ) were used for the tumor compartment and the side channel. The void fraction for modeling the delivery in the tumor-on-a-chip model was estimated by subtracting the volume occupied by CellTrace yellow-stained BT-474 cells in a Z-stack from top to bottom in the central tumor compartment from the total measured tumor compartment volume in the same Z-stack. The side channels are connected to the middle compartment by smaller channels that are 100 µm long, 20 µm × 20 µm wide, and 40 µm spaced. All of the boundaries were defined with a zero-flux condition, except for the side walls (normal to y-axis), which were set as symmetric. The initial conditions are shown in Figure 1B.
To determine the diffusion coefficient of DARPins in the ECM for the tumor-on-chip model, we performed simulations for the non-binding Off7 DARPin with different values for the diffusion coefficient. We then chose the value that provided the best fitting with experiments through a least-squares minimization method ( Figure S5A). The obtained value, of 55 µm 2 ·s −1 , was kept for the subsequent simulations with binding DARPins, assuming a comparable hydrodynamic radius for each.

Specificity of DARPins and Toxicity of DARPin-Toxin Fusion Proteins in 2D Model Systems
As a reference point for the evaluation in 3D tumor models, we started by investigating the binding specificity of DARPins and the cell-selective toxicity of DARPin-toxin fusions in flat 2D monolayer systems. We used DARPin 9_26, which binds to HER2 with a a dissociation constant (K d ) of 1.4 nM [35], and EpCAM-binding DARPin Ec1, which binds to EpCAM with a K d of 68 pM [36], including the non-binding control Off7 ( Figure S1A). DARPins were produced with a unique cysteine for fluorescent labeling or were fused to a fragment of ETA consisting of the translocation domain and the catalytic domain (Figure 2A and Figure S1B). DARPin-toxin fusions are highly potent receptor-targeted fusion toxins [39].
HER2-binding 9_26 and the EpCAM-binding DARPin Ec1 showed a clear membrane staining and some vesicular uptake when incubated with Calu-3 cells ( Figure 2B). The control DARPin Off7 showed very low non-specific signals. The specific binding of the DARPins was further quantified by flow cytometry (Figure 2C and Figure S3), which confirmed the highly specific binding of Ec1 and 9_26, but not of the control Off7. In co-cultures of IMR-90 fibroblasts and Calu-3 tumor cells, a specific binding of Ec1 and 9_26 to tumor cells was observed ( Figure 2D). No enrichment of the non-binding control DARPin on the tumor cells was seen.
Additionally, we evaluated the toxicity towards the cancer cell lines of DARPintoxin fusions containing either 9_26, Ec1, or Off7 after 48 h ( Figure 2E). We observed a 75% reduction in cell viability for both targeted fusion toxins, already at subnanomolar concentrations. To confirm that the residual viability was due to the short incubation period, we incubated cells for 48 h with fusions toxins, washed away the fusion toxins, and incubated them for an additional 72 h. At that time point, the targeted toxin mediated 100% tumor cell killing already at 10 pM, whereas the control fusion toxin showed no toxicity ( Figure S6).

DARPin Delivery in 3D Multicellular Tumor Spheroids (MCTS)
To evaluate the potential of mathematical modeling to predict DARPin diffusion in multicellular tumor spheroids (MCTS), we simulated DARPin delivery using a mathematical model similar to the one reported by Torres et al. for the investigation of in vivo protein delivery [28]. MCTS are a reliable 3D in vitro model for simulating the features of avascular tumors, including a binding-site barrier, because of a high receptor density. Optical clearing allows for the visualization of fluorescence throughout the entire spheroid using confocal microscopy [20,32]. Given our knowledge on most of the quantifiable parameters determining tumor penetration (Table 1), the unknown variable in our model for fitting variation in the speed of penetration is the apparent diffusion coefficient. We refer to this value as the apparent diffusion coefficient, because it takes into account the tortuosity of diffusion in the spheroid and is thus lower than the true diffusion. We simulated the tumor penetration of Ec1 and 9_26 at 500 nM and 100 nM in BT-474 spheroids with diffusion coefficients ranging between 0.5 to 10 µm 2 ·s −1 . For Ec1 at 500 nM, the simulations yielded a full penetration after~5.5 h with an estimated diffusion coefficient of 1.0 µm 2 ·s −1 (Table 1), which agreed with the experimental data, and at 100 nM, a full penetration of the spheroids would occur after~2 h at a diffusion coefficient of 10 µm 2 ·s −1 ( Figure 3B). For 9_26 at 500 nM, the model predicted that, contrary to Ec1, even with a diffusion coefficient of 2.0 µm 2 ·s −1 , complete penetration of the spheroids would take even more than 6 h ( Figure 3C). The different diffusion behaviors reflect mostly differences in receptor internalization rates. When simulating penetration at 100 nM for both DARPins, the binding-site barrier was more evident, showing that for 9_26, the simulated diffusion coefficients were not sufficient to reach the core of the spheroids within the 6 h time frame (Figure 3D), whereas for Ec1, only a diffusion coefficient of 10 µm 2 ·s −1 would achieve full penetration ( Figure 3B).
To experimentally validate our model, we investigated DARPin diffusion into tumors and observed similar time-dependent patterns of diffusion into the spheroids for both the EpCAM-targeting DARPin Ec1 and the HER2-targeting DARPin 9_26 at 500 nM ( Figure 3E,F). In a matter of hours, DARPins moved to the core of the spheroid, exhibiting the predicted "frozen front" at each time point, consistent with the presence of a bindingsite barrier [14]. After 6 h, Ec1 had fully penetrated the spheroid and 9_26 had reached about 90-95%. A lower concentration of 100 nM showed a slower-moving front, and not a generally lower signal ( Figure 3E,G). For Ec1, the simulations were most similar to experimental data assuming a diffusion coefficient of 1.0 µm 2 ·s −1 , whereas for 9_26, a simulation employing a diffusion coefficient of 2.0 µm 2 ·s −1 was most accurate for 500 nM, while at 100 nM, a faster apparent diffusion coefficient would be needed ( Figure 3A-D). 3D), whereas for Ec1, only a diffusion coefficient of 10 µ m 2 ·s −1 would achieve full penetration ( Figure 3B).  These results illustrate the potential for combining simulations of protein delivery in spheroids with optical clearing protocols for an accurate assessment of tumor penetration, allowing us to predict the effects of concentration and affinity on successful delivery, while also serving as an important validation step before moving to more complex microfluidic models.

DARPin Delivery in 3D in a Microfluidic Tumor-on-a-Chip Model
To be able to study DARPin delivery in real-time, we utilized a microfluidic tumor-ona-chip approach in conjunction with simulations. The ability to compartmentalize allows for precise control over the interface between the 3D tumor tissue and the drug solution, thus enabling an accurate real-time study of tumor penetration in a 3D tumor microenvironment. Here, we focused on transport by diffusion, mimicking tumor regions where little to no interstitial flow is present ( Figure 4A). In the microfluidic device, fluorescently labeled DARPins added in the medium reservoirs diffused from the side channels into the main tumor compartment containing BT-474 cells ( Figure 4B), and could be followed in real-time during 1 h using confocal microscopy ( Figure 4C). The modeling approach for the tumor-on-a-chip can be used to simulate the transport behavior and fit observed behavior in order to extract the apparent diffusion coefficients ( Figure 4D-G).
We estimated the diffusion coefficient for DARPins in the tumor tissue by fitting the diffusion profile of the non-binding DARPin Off7 to the simulated results of a range of diffusion coefficients ( Figure S5A), yielding an estimate of 55 ± 5 µm 2 s −1 , which is much higher than the diffusion constants in spheroids with high cell densities. In both the simulations and experimental data of Ec1 ( Figure 4D) and 9_26 ( Figure 4E), we observed slower transport into the tissue as a consequence of the predicted DARPinreceptor interactions. To assess the goodness of fit, we calculated the average relative error between the simulated and experimental curves for Off7, 9_26, and Ec1 ( Figure S5B-D). The fit was markedly better compared with a model without binding for 9_26 ( Figure S5E,F) and Ec1 ( Figure S5G,H). Evidently, there was no "frozen front" reflecting a strong binding-site barrier as observed in the spheroids. Instead, the transport behavior of the DARPins in the tumor tissue on-chip, with a lower cell density compared with spheroids, can be best described as hindered diffusion.
The model was further adapted for the delivery of agents that have toxic effects. Unlike fluorescently labeled agents that are constitutively endocytosed, the receptor regenerated cells that take up toxic agents will eventually be killed, thereby halting receptor regeneration and affecting the kinetics of the protein distribution. In the adapted model that assumed no receptor regeneration as a proxy for cell killing, we predicted the protein diffusion of 0.1, 1, and 10 nM of the targeted fusion toxins ( Figure 5A). Given the extremely potent toxic effects of DARPin-toxin fusions [39], we assumed that the remaining free receptor levels reflect the cell viability. The predictions indicated that extensive cell killing would occur only at 10 nM or higher. To compare the simulated data against the experimental data, we tested the activity of the DARPin-ETA fusion toxins in co-cultures of tumor cells with C5120 fibroblasts (1:1 ratio) in the microfluidic model. As our simulation predicted very little toxicity up to 1 nM, the co-cultures were incubated with 1 nM or 10 nM of Ec1-ETA DARPin and 10 nM of Off7 ETA-DARPin as a negative control. At 10 nM, but not at 1 nM, we observed a significant 10-fold increase in cell death, as reflected by the overlapping PI area in BT-474 cells, compared with the controls of Off7-ETA and the controls without DARPin ( Figure 5B), consistent with the prediction from our simulations. Furthermore, Ec1-ETA showed a specific killing towards high EpCAM-expressing BT-474 tumor cells, but not towards low EpCAM-expressing fibroblasts in the tumor-on-a-chip ( Figure 5B,C).

Discussion
In this study, we illustrated a hybrid approach that combines the simulation of protein behavior in 3D tumor tissues with an experimental investigation of tumor-targeted DARPins and DARPin-toxin fusion using spheroids and in real-time using a microfluidic tumor-on-a-chip model.
There is still a limited understanding of the dynamics of transport and binding interactions of protein-based targeted agents in complex tumor tissues, which is crucial for designing and evaluating novel therapeutic approaches for optimal target engagement in vivo [43]. We focused on spheroids that represent dense, solid avascular tumors with a high interstitial pressure, which often exhibit a layer-by-layer uptake of protein therapeutics and a binding-site barrier. Alternatively, we used tumors-on-chips to represent scenarios where the cells are more loosely associated and more prone to metastasize.
Mathematical modeling approaches have a long history in drug delivery research, though the large majority has been focused on nanoparticles [44,45]. In contrast, pharmacokinetic/pharmacodynamic modeling has long been used for predicting the overall receptor occupancy of protein-based therapeutics in order to guide dose selection in humans, but these provide little mechanistic insights at the level of the tumor tissue [46].

Discussion
In this study, we illustrated a hybrid approach that combines the simulation of protein behavior in 3D tumor tissues with an experimental investigation of tumor-targeted DARPins and DARPin-toxin fusion using spheroids and in real-time using a microfluidic tumor-on-a-chip model.
There is still a limited understanding of the dynamics of transport and binding interactions of protein-based targeted agents in complex tumor tissues, which is crucial for designing and evaluating novel therapeutic approaches for optimal target engagement in vivo [43]. We focused on spheroids that represent dense, solid avascular tumors with a high interstitial pressure, which often exhibit a layer-by-layer uptake of protein therapeutics and a binding-site barrier. Alternatively, we used tumors-on-chips to represent scenarios where the cells are more loosely associated and more prone to metastasize.
Mathematical modeling approaches have a long history in drug delivery research, though the large majority has been focused on nanoparticles [44,45]. In contrast, pharmacokinetic/pharmacodynamic modeling has long been used for predicting the overall receptor occupancy of protein-based therapeutics in order to guide dose selection in humans, but these provide little mechanistic insights at the level of the tumor tissue [46]. Very recently, Tang et al. reported an investigation of the dynamics of antibody-target binding in tumor tissue in vivo. Notably, they reported large differences between k on and k off , as measured by surface plasmon resonance (SPR) and apparent rates obtained after fitting the model with in vivo data, as well as very large differences in apparent k off rates between stroma-rich and stroma-poor tumor regions. As the model was compartmental, diffusional hindrances were not explicitly modeled in their approach, which could have contributed to the observed slow on-and slow off-rates in the stroma-rich tumor region, as also indicated by the authors [47]. Their approach differs from our approach, where we assumed fixed k off rates and fitted the data according to apparent diffusion coefficients in a spatiotemporal fashion, allowing for an analysis of the gradients.
Our modeling approach is more similar to the approaches followed previously by Thurber and Wittrup, who performed a spatiotemporal investigation of the tumor delivery and binding-site barrier of targeted proteins in vitro [14], and Torres et.al, who performed a theoretical investigation of in vivo protein delivery [28].
In agreement with their findings, we observed a binding-site barrier on the outside of the spheroids when using tumor cells that overexpress membrane molecules such as EpCAM or HER2, and DARPins against these targets. The "frozen front" moves in predictable ways as a function of concentration and time, as has also recently been found in other contexts by us and several other groups who have studied the penetration of binding proteins in spheroids [14,32,48,49]. When fitting the experimental data of the DARPin diffusion in the spheroids, the model predicted diffusion coefficients of 1 µm 2 ·s −1 for Ec1 and 2 µm 2 ·s −1 for 9_26. It must be noted that these values rely on reported internalization rates in 2D cultures, which may be different from those in 3D tumor tissues [37,38]. It is noteworthy that the diffusion coefficients we found were low when compared with the 33 µm 2 ·s −1 reported by Thurber and Wittrup for a single-chain antibody fragment in tumor tissue [14]. Factors that could decrease the diffusion coefficient of the DARPins in the spheroids are related to the density of the collagen or methylcellulose within the spheroids, the cell density and the tightness of the cell-cell junctions [50]. In Thurber and Wittrup's system, the same spheroid size was obtained with ten times fewer cells, which may have yielded more loosely associated spheroids in which proteins can diffuse faster. In the microfluidic tumor-on-a-chip model, tumor cells grow at a much lower density, explaining the much higher estimated diffusion coefficient of 55 ± 5 µm 2 ·s −1 . This latest value is similar to what would be expected in a relatively dense tumor tissue in vivo for a similarly sized molecule, reflecting, for instance, U87 and HSTS26T tumors grown in vivo in a mouse dorsal chamber [51]. Of note, even though the Matrigel used in this study is mostly a mixture of laminin, collagen IV, and enactin, and reflects a basement membrane [52], while in tumor tissues it is mostly collagen that lowers the diffusion coefficient [51], the transport properties of DARPins in cell-dense tumor tissues can be investigated properly when diffusion coefficients match what is found in vivo and specific binding to matrix components plays a negligible role, as is expected for hydrophilic DARPins with highly specific recognition surfaces.
In contrast with the spheroids model, the tumor-on-a-chip system is a model that allows for the investigation of tumors with lower cell densities, in which the accessibility depends, to a large degree, on the porosity of the ECM that causes diffusional hindrance [51]. In these kinds of tumors, the proteins will face diffusion through the extracellular matrix as a main challenge, a notorious problem especially in pancreatic adenocarcinomas [53]. Hence, the proteins will experience a hindered diffusion through the 3D matrix, which is in line with the simulations of our mathematical model, and which is distinct from the frozen front observed by us and others in spheroids. The developed model was also employed to simulate the transport of a DARPin fused to a toxin, which specifically binds and selectively kills tumor cells in the microfluidic device. This combined approach can lead us to better study the specific delivery of protein therapies in a complex 3D system, while also predicting the minimally needed concentration for optimal therapeutic effects. Further applications of our approach can bring us closer towards rational therapeutic protein design. Combined with developments in high-throughput production and approaches towards de novo protein design with desired properties [54], faster in vitro testing using microfluidic organ-on-chip technologies could help with predicting and developing novel and more powerful therapies that can be tested more rapidly and be faster translated to the patient [55].

Conclusions
To summarize, the key achievement of this study is the illustration of the power of a hybrid approach that combines the use of biomimetic in vitro tumor models with mathematical models to study and predict protein delivery and activity in complex 3D systems. Our hybrid approach has great potential for the development of novel potential protein therapies, as it in part sidesteps the extrapolation challenges from 2D or animal models to complex human physiology by efficiently combining in silico approaches with biologically relevant in vitro test systems.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cancers13102461/s1, Figure S1: Purified DARPins and DARPin-toxin fusions, Figure S2: Flow cytometry-based quantification of receptor density in different cell lines, Figure S3: Binding of DARPins to BT-474 cells, Figure S4: Schematic representation of the input variables and outcomes determined by the binding-diffusion mathematical model, Figure S5: Fitting of the mathematical model to the experimental data, Figure S6: Long term viability assay after incubation with DARPintoxin fusions, Supplementary methods, 1: quantification of cell death in tumor-on-a-chip system, 2: analysis of spheroid penetration. Funding: This research received the support of Radboudumc, through an internal funding program, and the support of Institut Pierre-Gilles de Gennes (laboratoire d'excellence, "Investissements d'avenir" program ANR-10-IDEX-0001-02 PSL and ANR-10-LABX-31).

Data Availability Statement:
The data that support the findings of this study are contained within the article and supplementary materials. Raw data are available from the corresponding author upon reasonable request.