The tumor stroma is the abetting structure of tumor tissue. It is composed of non-malignant cells of tumor-like cancer-associated fibroblasts (CAFs) (Raghu et al., 2016; Tang et al., 2016; Lee et al., 2018; Liu et al., 2019; Monteran and Erez, 2019; Wang et al., 2019; Yavuz et al., 2019), tumor-associated macrophages (TAMs) (Qian et al., 2009; Noy and Pollard, 2014; Larionova et al., 2019; Lin et al., 2019), tumor-associated neutrophils (TANs), mesenchymal stem cells and extracellular matrix (ECM) with fibrous structural proteins (e.g., collagen and elastin), fibrous adhesive proteins (e.g., laminin and fibronectin), and proteoglycans. Nests of malignant tumor cells are linked through junctional proteins (e.g., claudins, desmoglein-2, and E-cadherin) in most solid tumors derived from epithelial tissues. These nests are surrounded by tumor stroma which plays a key role in regulating the behavior of cells found in the local milieu. Tumor stroma creates a niche that aids in seeding metastatic cells and intervenes in drug delivery to tumors. It generates a physical obstacle of stroma proteins that restrict drug penetration and connection between drug candidates, tumor-infiltrating immune cells, and target receptors present in the tumor cell surface. Tumor stroma generates cytokines and chemokines, which prompt synthesis of stroma proteins, prevent activation of immune cells, and activate immuno-suppressive cells such as regulatory T cells. Stroma is associated with characteristic tissue development and homeostasis in the tumor microenvironment, distinct from those associated with normal tissue (Kim et al., 2013; Nunes et al., 2019). Moreover, ECMs produced in most tumors make them more rigid and different types of collagen molecules are also accumulated, forming a thick network inside tumor tissue, resulting in a decrease in pores of tumor matrix, which restrict tumor penetration of nanomedicine useful for therapy. Amplified rigidity increases interstitial fluid pressure, which further restricts the distribution of nanomedicine throughout tumor mass. Therefore, targeting genetically stable stromal cells provides an additional benefit.

Stromal cells exhibit constant synthesis and release of growth factors, connective tissue components, and cytokines, cooperate with malignant cells to proliferate, invade, and metastasize, which seek major attention in tumor-targeted drug delivery using nanomedicine (Choi et al., 2013; Raghu et al., 2016; Tang et al., 2016; Lee et al., 2018; Liu et al., 2019; Monteran and Erez, 2019; Wang et al., 2019; Yavuz et al., 2019). Cancer-associated fibroblasts (CAF) are major cells found in tumor-associated stroma compared to stroma cells in healthy tissue (Kalluri, 2016). CAF are spindle-shaped mesenchymal cells characterized by constant activation, faster proliferation, and accumulation of large amounts of ECM compared to fibroblast in normal tissue (Kim and Bae, 2016). Tumor cells and stromal cells upregulate different types of profibrotic growth factors in the tumor microenvironment to transdifferentiate stromal fibroblast in CAF. They release various growth factors, for instance, epidermal growth factor (EGF), hepatocyte growth factor (HGF), and insulin-like growth factor-1 (IGF-1) and affect cell proliferation, invasion, and metastasis. They are involved in angiogenesis along with inflammatory cell recruitment. CAF are also engaged in the arousal of angiogenic processes and engages more inflammatory cells. Another very important cell type present in tumor stroma is tumor-associated macrophages (TAMs), which are immune cells of the tumor microenvironment. TAMs suppress antitumor immune responses, generate an immune suppressive microenvironment, and control tumor-associated changes in ECM by secretion of profibrotic growth factors. They produce cytokines (IL-I and IL-8), tumor necrosis factor-α (TNF-α), growth factors (EGF, HGF, bFGF, and VEGF), and various enzymes (Li et al., 2017a). TAMs regulate cancer stem cell activities in solid tumors. Tumor-associated neutrophils (TANs) are also a dominant form of immune cell infiltrates, found in various types of cancer (Masucci et al., 2019; Wu et al., 2019). They generate neutrophils, reactive oxygen species, cytokines, growth factors, and proteinases and play key roles in controlling tumor cell proliferation, metastasis, angiogenesis, and antitumor immune suppression (Masucci et al., 2019; Wu et al., 2019).

Multicellular tumor spheroids are spherical in vitro self-assembly of cellular aggregation representing the characteristics of in vivo solid tumors. In this self-assembled organization, cells aggregate, sort, and compartmentalize to separate different regions of the spheroid (Nath and Devi 2016). To impersonate the shear properties of the in vivo model of solid tumors, the study of the internal structure of the in vitro model of a spheroid is essential. The internal structure is incorporated with different cell layers based on the concentration gradient of nutrients, oxygen to regulate cell function, differential zones of proliferation and growth factors, pH, and cellular density. All these factors of the cellular organization play an indispensable role in studying the diffusion and penetration for impaired nanomedicine delivery (Lazzari et al., 2017; Sant and Johnston, 2017b; Costa et al., 2016). In accordance with differential proliferation rate, the structure can be broadly categorized in three different zones: proliferation, senescent, and necrotic zones; see Figure 5 (Mehta et al., 2012; Lazzari et al., 2017; Sant and Johnston, 2017). Modeling of spheroid growth can be done by taking into account the following spheroid cellular organization. These layers are characterized by the decreasing gradient of nutrients, oxygen, and pH, from the exterior to the center of the spheroid and by the increasing gradient of CO₂, lactate, and waste from the exterior to the center of the spheroid. The outer layer where cells resurge and escalate rapidly in number is called the proliferation zone. It exhibits high proliferation rates in the spheroid periphery. The proliferation is stimulated by the constant exposure of the cells to oxygen and nutrients. When the oxygen diffusion and nutrient availability become a limiting factor, the cell’s proliferation rate decreases, giving rise to the middle layer made up of senescent cells. In the middle zone, cells can no longer divide due to the depletion of nutrients; nevertheless, they are active and alive, which is identified as the senescent zone. The middle region is followed by a decrease in cell metabolism as the distance from the outer region increases. The supply of nutrition is also depleted as we move from the periphery to inner zones. Finally, the core region, characterized by depletion in oxygen concentration, results in hypoxia, nutrient supply, and waste accumulation leading to a critical situation of cell necrosis (Figure 5). It is identified as the region where cells are noxious and lead to death, known as the necrotic zone. This core is the zone of the lowest pH (6.6–7.2) within the tumor spheroid. In this hypoxic environment, the pyruvate is converted into lactate by the cancer cells to obtain energy, known as the Warburg effect. The accumulation of lactate results in an acidified core of the spheroid and makes it favorable for drug release from nanomedicine.

The annular organization of subsequent cell layers within the spheroid is reported to be the key factor behind the inadequate therapeutic efficacy of anticancer drugs, as the presence of cells in the successive region from the outer zone, that is, the quiescent and necrotic region, reduces the potency to penetrate nanomedicine. In the hypoxic region of the spheroid, different drugs that cause cell death via the formation of reactive oxygen species (ROS) have a very low therapeutic efficacy (Liang and Liu, 2016; Al-Akra et al., 2019). Moreover, many drugs (e.g., paclitaxel) that promote cytotoxicity to the cells and proliferation have a poor impact in this innermost zone and that is why even the drugs which are highly effective in rapidly multiplying cells have a poor therapeutic benefit in the interior layers of the spheroid for the same reason. Moreover, low environmental pH catalyzes the change of net charge of drug molecules (e.g., melphalan, methotrexate, and mitoxantrone), tumor penetration, and intracellular uptake. The typical low pH condition at the interior zone induces the cancer cells to start expressing different kinds of hypoxia-inducible factor (HIF) family, which promotes the proliferation and consequent survival rate of the cancer cells (Masoud and Li, 2015; Albadari et al., 2019; Hompland et al., 2021).

The application of multicellular tumor spheroid in tumor-drug delivery is increasing rapidly. The key characteristics of a promising in vitro spheroid model are to accurately mimic the in vivo solid tumor, incorporating the real biological characteristics such as heterogeneity, microenvironment, and cellular arrangements. Spheroid is a well felicitous in vitro model for the aimed conveyance of nanomedicine due to its ability to transcribe the intrinsic characteristics found in solid tumors. We will elucidate the following characteristics of the spheroid in detail such as cellular heterogeneity, cell-cell signaling, internal structure, ECM deposition, gene expression, and how the external modification to spheroid could lead to efficient drug delivery, and how mathematical modeling approaches could give us the tools to decipher the underlying mechanisms (Albrecht et al., 2006; Ekert et al., 2014; Popilski and Stepensky, 2015; Shamsi et al., 2019).

Three-dimensional multicellular tumor spheroid is widely used as a model system for studying different types of cancers and screening of nanomedicine efficacy. Cellular heterogeneity is one of the critical aspects that is needed for this study and represents the diverse nature of cell types showing distinct morphological and phenotypic profiles. Cellular heterogeneity is embodied with the multiple cells that vary in their protein and gene expression. Spheroid is embraced with different heterogeneous tumor and stromal cells with varying ratios that lead to its cellular heterogeneity. The in vitro model can be established with homotypic cancer cells only or heterotypic cancer cells cultured with fibroblast, immune cells, or endothelial cells (Costa et al., 2016; Lazzari et al., 2017; Millard et al., 2017). In heterotypic tumor spheroids, the number of cells could be varied to closely mimic heterogeneity in solid tumors. Moreover, interacting networks among cells having this heterogeneity are demanded. Signaling established between the cells which are in close proximity helps in reproducing the spheroid’s key aspects and its microenvironment. All the cells growing in close proximity provide a route to a highly interactive microenvironment for physical communication and signaling pathways of cells to drugs. This heterogeneity and signaling can help decipher how multiple cell types might impact drug delivery to the site.

Solid tumors exhibit complex cellular microenvironment and architecture, where extracellular matrix (ECM) plays a crucial role (Lu et al., 2012; Pickup et al., 2014; Bussard et al., 2016; Lazzari et al., 2017). ECM is an essential component for binding the cells together via proteoglycans and fibrous proteins. The extracellular matrix (ECM) compositions, e.g., fibronectin, laminin, proteoglycans, collagen, and tenascin, are deposited by cells within 3D spheroids like solid tumors (Frantz et al., 2010; Kim et al., 2011; Lu et al., 2012; Hynes, 2013). The ECM (α5- and β1-integrin) in the tumor spheroid form a barrier known as limited mass transport, further limiting the insertion and uniform distribution of nanomedicine in the solid tumor mass. Extracellular matrix components in a highly interactive 3D microenvironment of spheroids make the structure more compact by depositing the protein, leading to an increase in the interstitial fluid pressure and further limiting impaired therapeutic efficacy of nanomedicine. It helps tumor cells in regulating different cell functions and maintaining the complex cell network. ECM-related signaling pathways play a crucial role in tumor progression, cancer cell migration, and penetration (Frantz et al., 2010; Kim et al., 2011; Lu et al., 2012; Hynes, 2013). Tumor cells within solid tumors grow in two stages. Initially, the tumor grows very rapidly, termed the avascular growth phase. Then, cells continue to stay in a dormant phase followed by vascular growth with the generation of new blood capillaries (angiogenesis) promoted by pro-angiogenic factors (e.g., angiogenin) and ECM rebuilding mediated by proteases (e.g., MMP-2 and MMP-9). Similarly, the progress of a tumor in an in vitro model is the collective dynamics of interaction between tumor cells and their microenvironment. The initial volume growth increases exponentially with time until it reaches a certain value (∼400 μm in diameter), but later the growth of volume decreases with time and becomes constant, referred to as the spheroid growth plateau (Figure 5). This avascular tumor expansion, growth plateau, and vascular expansion as spheroid growth varying with both temporal and spatial dimensions can be modeled mathematically.

Growth factors and the cellular ECM protein can be encoded with the learning of genes. Gene expression is a closely constrained process regulating the response of a cell to its changing environment. The study of gene expression is vital as it involves the conversion of DNA’s instruction into functional products like protein (Costa et al., 2016). Gene expression profile is firmly affected by the cellular organization of the spheroid. The abundant target protein expression, also known as overexpressed gene, plays an essential role in studying cancer progression, invasion, and metastasis. There are various factors that influence the morphology of tumor spheroids, such as cell type, cell density, culture media, method of culture, and mechanical stress. Based on their compactness, spheroids could be compact spheroid, tight aggregates, and loose aggregates of cells. In a compact spheroid, cells are tightly bound to each other. In tight or loose aggregation, cells do not form a complete sphere and can be easily disintegrated. Aggregation of cells initially occurs by loose bonds between integrin and ECM, then forming close contact through N-cadherin to E-cadherin interactions. As the cell communication pathways, morphology, and polarity of the cells in solid tumors closely resemble the structure of multicellular tumor spheroids in many ways, the study of using tumor spheroids is a very important in vitro model in recent times for studying the penetration profile of nanomedicine and calculation of accumulation of anticancer drugs (England et al., 2013). Thus, cellular arrangements and internal structure in tumor spheroid closely resemble solid tumors, making it an appropriate model for the study of tumor growth and invasion as metastasis with the effect of drug candidates for screening of different nanomedicine formulations towards its efficacy and efficiency of disease management.

A plethora of techniques to culture in vitro models are available, incorporating the use of cell attachment resistant surface forces to induce cell-cell interactions and ultimately support the formation of multicellular tumor spheroid (Figure 3 and Table 2). These techniques will be discussed in the subsequent part (Ishiguro et al., 2017; Lazzari et al., 2017).

Advanced characterization techniques are utilized to characterize tumor spheroids based on, e.g., size, shape, cellular arrangements, protein and gene expression, cell cycle patterns, invasive nature, and metastatic potential of cancer cells to assess the nanomedicine (Figure 3 and Table 2). Various types of techniques are described as follows (Elizabete et al., 2016).

The study of spheroid growth using mathematical modeling is more than half a century old (Laird, 1964). These quantitative approaches are invaluable tools for comprehending the cellular and transport phenomenon within spheroids and foreseeing the physiological acknowledgment to drug delivery. They can provide the pertinent perspective for drug delivery mechanisms and effectiveness in spheroid. Mathematical modeling of spheroid can be outlined into two distinct scales. First, the cellular scale describes cellular dynamics that lead to the model of tumor cell proliferation. The exact quantitative expression for this model is given by the Gompertz equation (Laird, 1964). It was commonly believed that tumors grow exponentially and stop until the host nutrition supply is exhausted. However, it has been observed that exponential growth is only dominant for a brief period of time and reaches a growth plateau as the tumor grows larger in size (Laird, 1964; Ward and King, 1999; Mehta et al., 2012; Costa et al., 2016). Secondly, on the other hand, the subcellular and extracellular scales describe the chemical dynamics with the help of the reaction-diffusion mechanism and hybrid models. Modeling tumor morphology starting from individual cells is usually divided into two categories incorporating both continuum and the cell-level description. Hybrid models deal with a combination of these two different descriptions. One is related to the periphery of the tumor, embraced with a cell-level description where it is advantageous to do so, and the other covers the two inner zones of tumor and the extracellular matrix pertaining to continuum description, i.e., cell population-based continuum models and individual cell-based discrete models. During the development of these models, the primary focus was on incorporating the different characteristics of spheroid growth, such as initial exponential growth and the concentration gradient of nutrients, oxygen, which are vitally important for the layered organization of spheroid and targeted nanomedicine delivery.

First, cellular dynamics can be studied mathematically by one of the pioneering models for tumor growth developed by Laird in 1964 (Laird, 1964). The growth kinetics of solid tumors is akin to spheroid, which can be classified into two levels. During the initial phase, the exponential growth of the tumor volume is observed. It is followed by a dormant phase of minimized metabolic activities until the spheroid growth plateau, where the spheroid’s volume attains a constant value (Figure 5). In improving the exponential model, the cell population growth curve with a time-dependent growth rate is considered. Let the size of the population at time t be W ( t ) and the growth rate decay exponentially be a ( t ) = α e − b t . Here, the independent variable is time t , and the dependent variable is the tumor volume or population size W ( t ) . The corresponding ordinary differential equation for W ( t ) can be written as follows: d W ( t ) d t = α e − b t W ( t ) . (1)

The solution to this model shows tumor cell proliferation that can be expressed by a modified exponential process, commonly known as the Gompertz equation for sigmoidal growth, of the following form: W ( t ) W 0 = e a b ( 1 − e − b t ) , (2) where W ( t ) is the tumor size at any time t , W 0 is the initial tumor size, and   b is constant. Now, e − b t can be expressed in power series as   e − b t = ∑ n = 0 ∞ ( b t ) n n ! . During the initial growth, i.e., at small t , e − b t can be approximated as e − b t ≈ 1 − b t , and the growth equation takes the simple exponential form W ( t ) W 0 = e a t , which is consistent with the observation of initial tumor growth. At later times, the growth deviates from the pure exponential dependency and takes the Gompertz form. From the Gompertzian analysis, the theoretical upper limit of tumor growth for mice is typically ≈ 10 9 cells, which is also consistent with the approximate size at death (Laird, 1964; Norris et al., 2006; Altrock et al., 2015; Costa et al., 2016).

In the above formalism of growth dynamics, there is only one independent variable time t , and the dependent variable is the volume of the tumor. Nevertheless, mathematical models have been developed using partial differential equations to study the spheroid growth and architecture at higher dimensions. For instance, a model has been developed to study quantities such as oxygen distribution which has more than one dependence, one is spatial and the other is temporal. One valuable characteristic of a spheroid is their limit of diffusion of about 150–200 μm for many molecules, specifically oxygen (Mehta et al., 2012). This diffusion limit gives rise to limited mass transport, as a result of which spheroid displays the gradient of oxygen, distribution of nutrients, metabolic waste accumulation, and proliferation profile inside them. Hence, a diffusion model for the oxygen concentration is requisite. The oxygen concentration in a tissue can be represented by ρ ( r , t ) , where the position r , in general, measures from the center of the spheroid at a time t. The reaction-diffusion equation for oxygen concentration inside a spheroid in the form of a partial differential equation can be written as follows: ∂ ρ ( r ,   t ) ∂ t = D ∂ 2 ρ ( r ,   t ) ∂ r 2 − a ρ ( r , t ) − b N ( r , t ) ρ ( r , t ) + c ρ ( r , t )   , (3) where oxygen diffuses inside the spheroid with a diffusion constant D , decays with the rate constant a , and is produced with a rate constant c . The consumption of oxygen is proportional to the group size of the tumor cells N ( r , t ) and available oxygen concentration ρ ( r , t ) itself with a rate constant b .

A mathematical model is required to weave the insights gained from the discrete cellular dynamics into a coherent description of the reaction-diffusion mechanism. Combining the continuum model of growth-consumption as a reaction-diffusion model along with discrete cellular dynamics of cell growth and motility, a hybrid model has been developed. The model retains the cellular description in the rapid proliferation region on the periphery of the tumor and for the dynamics of tumor cell density, extracellular matrix (ECM) cells, matrix-degrading enzymes (MDE), and oxygen concentration as continua. The corresponding coupled dynamic equations of these individual quantities can be written as partial differential equations. For instance, a partial differential equation for the dynamics of tumor cell number N T is given by ∂ N T ∂ t = D T ∂ 2 N T ∂ r 2 − ∂ ∂ r ( N T ∂ N E ∂ r ) , (4) where the diffusion constant of tumor cells is denoted by D T   and extracellular matrix is represented by N E . The dynamics of ECM ( N E ) , MDE ( N M ) , and oxygen concentration ( ρ ) are represented by ∂ N E ∂ t = − ϕ   N M N E , ∂ N M ∂ t = D M ∂ 2 N M ∂ r 2 + λ   N T − μ   N M , ∂ ρ ∂ t = D ρ ∂ 2 ρ ∂ r 2 + f   N E − n   N T − c ρ   . (5)

The diffusion constant associated with ECM and MDE is, respectively, represented by D E and D M . The degradation of the extracellular matrix is directly proportional to the density of matrix-degrading enzymes and extracellular matrix with a proportionality constant ϕ . Matrix-degrading enzymes are produced by the tumor cells themselves with a rate constant λ and natural decay with a rate constant μ . On the other hand, the oxygen concentration is directly proportional to the ECM density with a rate constant f . It is consumed by the tumor cell at a rate n and decays at a rate c .

The hybrid model bridges two different types of models distinguishing between individual cell and cell population-based models. This model complements a fully continuous and complex description of tumor dynamics. The discrete cell interaction can be explained by a stochastic model, which is an off-lattice model. On a 2D lattice, probabilistic rules are applied to each cell pertaining to discreteness by defining stochastic reaction rates of respective events and may depend on the microenvironment (Kim and Stolarska, 2007; Altrock et al., 2015). In a hybrid model, these stochastic rates depend on the concentration of continuous variables ρ , N M , N E , and N T   .   Besides this, the hybrid model also specifies the guidance for processes like proliferation, which are dependent on the environment and are specific to cells. Recent research has been done considering 3D cell culture leading to new methods for drug transport. In the future, mathematical models will continue to help as a guiding path for understanding the tumor architecture and its growth and studying the transport of oxygen and nutrients among different zones within the spheroid via cellular and chemical dynamics.

The organization of interior layered composition of MCTS can be determined by the spatial arrangements, interaction, and grouping of cells combined into developmental and functional patterns. The fabrication of these similar or different cellular patterns allows them to form self-organized individual compartments assembled together, resulting in a highly stratified structure is known as self-assembly. This bio-fabricated spheroid assembly is mediated by the molecular gradients of soluble factors which are capable of binding to cellular receptors causing the signals to initiate proliferation. These soluble factors within the spheroid microenvironments are established by the process of convection and diffusion. This arrangement across the multiple length scales can be assembled in two different ways: cluster-based self-assembly and collision-based self-assembly. Cluster-based self-assembly consists of the partitioning of mono-dispersed cells into sectors followed by their settlement and aggregation, maintaining their cellular integrity as clusters and resulting in spheroid. In contrast, when the suspended cells strike into each other, leading to the formation of a spheroid is termed collision-based arrangement. Self-assembly is often followed by self-sorting. When varying cell types are organized among themselves, leading to a particular pattern of segregation, this can be stated as self-sorting (Figure 6). Based on this, a theoretical model has been developed for the study of spheroid formation that occurs when two cell types are differently segregated. Modeling of a spheroid based on self-assembly and sorting is affiliated with cell-to-cell adhesion and surface tension. This necessitates the sorting of cells of the highest cohesion to the interior of the spheroid and those with lower cohesion to the outer. This hypothesis is known as Differential Adhesion. From the recent study (Achilli et al., 2012), it has been shown that modeling of self-sorting processes can be done using an order parameter to describe the relationship between heterotypic interface length and size of the system, taking the geometrically driven argument as its base.

The subsequent part of the study will elucidate another key component of the targeted drug delivery mechanism, that is, the study of drug molecules and how they can be modeled mathematically and the key attributes related to their physical properties and chemical composition.

The rapidly growing nanotechnology evokes the need for understanding how nanomedicine’s characteristics influence the transport processes. Mathematical modeling has the potential to point to the comprehensive view of nanomedicine designing and characterize the important features prerequisite for drug delivery. In this section, we will discuss the quantitative insights for the elucidation of nanoparticle (i.e., nanomedicine) diffusion mechanisms in the bulk, penetration into the multicellular spheroid, and then the calculation of binding site availability on the spheroid surface (Figure 7). The viability of requisite mathematical models for the treatment of spheroids with nanomedicine will be of primary focus.

The binding of a drug molecule is essential when it comes to an unhindered drug delivery mechanism, which is proportional to the optimized value of successful binding. Molar concentration plays a major role while addressing the optimized binding and penetration of nanoparticles to and through the spheroid. Binding can be maximized in two different ways: first, it is related to the molar concentration of nanoparticles, while the other is related to the molar concentration of the cell surface binding sites in the spheroid. To study the penetration of nanomedicine to the spheroid, using stock’s solution (Gao et al., 2013), the molar concentration of nanoparticle ( ρ M ) in Mol/L can be calculated as follows: ρ M   ( M o l / L ) = N P   c o n c .   ( n u m b e r / L ) 6 × 10 23 p a r t i c l e s / m o l e = 6 C s t o c k   × 10 12 N A   π ρ d 3 , (6) where C_stock is the nanoparticle concentration in manufacturer stock solution, N A is Avogadro’s number, ⍴ is the density of nanoparticle, and d is the nanoparticle diameter in meters. The stock’s solution is being prepared for the transport study of the molar concentration and aggregation of nanoparticles. Moreover, the molar concentration of cell surface binding sites in a spheroid can be defined as follows: B m a x ,   s p h e r o i d = ϕ × 10 15 6 × 10 23 × 4 3 π r c e l l 3 B m a x ,   s i n g l e   c e l l . , (7) where ɸ is the cell density calculated by taking the volume of total cells in spheroid divided by the total spheroid volume (with the assumption of the spherical shape of spheroid), r_cell is the radius of the tumor cell in the suspension of stock’s solution, B_{max, single cell} is the molar concentration of single-cell binding, and B_{max, spheroid} is the molar concentration of cell surface binding site in spheroid (Gao et al., 2013).

Mathematical modeling is instrumental in aiding our understanding of the internalization of nanomedicine through the stratified structure of tumor spheroid and then followed by its diffusion inside the spheroid. The diffusion of nanomedicine at the cellular and extracellular levels is stipulated by its calculation in the interstitium and interstitial sites and taking the porosity of the spheroid into account, addressed using mathematical tools. Nanoparticle penetration to the cells can be designated as a two-step process. First, the adsorption of nanoparticles is then followed by their internalization. Fundamentally, association or this binding kinetics of negatively charged nanoparticles to cell surface termed Langmuir adsorption with the assumption that a particular cell line has a fixed capacity for particle binding through electrostatic interactions (Wilhelm et al., 2002). This process incorporates the rate of change of mass of nanoparticles bound to the cell surface, established by considering the molar concentration in the extracellular medium (M), the mass of adsorbing and desorbing nanoparticles. Langmuir adsorption can be modeled using a differential equation representing the variation of m(t) with time consisting of k a s s , k d i s s   as the association and dissociation rate constants: d m d t = k a s s M ( m 0 − m ) − k d i s s m , (8) where m 0 and m ( t ) are the maximum mass that can be bound to the cell surface initially and at a later time t, respectively. The above differential equation represents that the mass of absorbed nanoparticles per unit time is proportional to the molar concentration of nanoparticles, to the mass that is present on the cell surface that still can be adsorbed ( m 0 − m ) , and to the mass of desorbing nanoparticles. Initially, at t = 0, m(t) = 0, the analytical description of mass of the adsorbed particles unifying both temporal dependencies and the concentration of binding nanoparticles on the cell surface can be given by m ( t , M ) = k a s s M k a s s   M +   k d i s s m 0 ( 1 − e − ( k a s s M + k d i s s ) t ) . (9) The above expression can be modified to get maximum adsorbed mass with characteristic time   T = 1 ( k a s s M +   k d i s s   ) , which approximately gives. m m a x = k a s s M k a s s M + k d i s s m 0 .

After adsorption has occurred, the internalization of nanoparticles can be quantitatively approached in two ways. First, by considering the mass of the nanoparticle bound to the cell surface and another by taking the structural uniformity of the target, i.e., tumor spheroid, into account. Starting with the first approach, the global process of nanoparticle penetration into the tumor involves their binding mass to the cell surface represented at any time t by m b i n d ( t )   and their internalized mass within the cell via endocytosis represented by m i n t ( t ) . Differential equation regulating all these phenomena can be stated as follows: d m b i n d ( t ,   M ) d t   = k a s s M ( m 0 − m b i n d ) − k d i s s m b i n d − d m i n t d t . (10)

It has been previously assumed that the fraction of cell surface absorbs the nanoparticle, termed reactive surfaces (RS) which are regenerative and remain constant. Let f i n t ( t ) ,   f 0 be defined as the fraction of RS being internalized over the total available reactive surface at any time t and the maximum fraction of RS internalized by the cell, respectively. Thus, this can be represented by the differential equation with internalization rate constant ( k i n t ) , d f i n t ( t ) d t = k i n t ( f 0   − f i n t ) . This equation can be used to calculate the rate of change of internalized mass as d m i n t d t = d f i n t d t m b i n d . Using this equation, we can also calculate the total uptake of the mass of the nanoparticle m ( t ) = m b i n d + m i n t .

Now, considering the structural uniformity, the penetration of nanomedicine into the spheroid is addressed quantitatively. Mathematical models have been developed using the nanoparticle-cell bio interface data to predict the penetration of nanomedicine into the spheroid that accounts for radial dependent changes of its internal structure. The diffusive transport and the spatial distribution of nanoparticles can be represented using nanoparticle-cell interaction parameters: association, dissociation, and internalization rate constants. The mathematical expression for the nanoparticle motion into the spheroid relationship is given by Goodman et al. (Goodman et al., 2008). The kinetics of nanoparticle internalization in 3D multicellular spheroid depends on the molar concentration of free nanoparticles available in the spheroid, bound nanoparticles, nanoparticles undergoing internalization, and concentration of available binding sites on the cell surface: M, M_b, M_int, and M_bs, respectively. The entire ensemble of defining free, bound, and unbound nanoparticles is dependent on the force, which in turn is potentially mediated by cell-cell interaction and cell-ECM interaction. A set of partial differential equations representing the rate of change of respective concentration with time taking the structural changes of the spheroid in the radial direction and nanoparticles internalization using concentration into account is given (Graff and Wittrup, 2003; Ward and King, 2003; Hori et al., 2021; Gao et al. 2013). The molar concentration of free particles per unit time in spheroid volume is proportional to the molar concentration of unbound particles present in the accessible volume intracellularly, which are varying with radial coordinates, dissociated bound particles, and associated binding sites as well as unbound particles concentration can be written as follows: ∂ M ∂ t     = D r 2 ∂ ∂ r [ ϵ r 2   ∂ ∂ r ( M ϵ ) ]   + k d i s s M b − k a s s M b s M ϵ , (11) where k_ass, k_diss, and k_int are the kinetic reaction rates for nanoparticle-cell interaction corresponding to the association, dissociation, and internalization of nanoparticles in a spheroid, respectively. D is the effective diffusion coefficient, r is the radial coordinate representing nanoparticle diffusion measured from the center of the spheroid, R is the spheroid radius, and ϵ is the volumetric porosity of the spheroid, which is the ratio of the spheroid volume accessible to the particles to the total available volume. The partial differential equation elucidating the change of molar concentration of bound nanoparticles with time embedded proportionality with association rate of binding sites, dissociation, and internalization rates of bound particles can be presented as follows: ∂ M b ∂ t   = k a s s M   b s M ϵ − ( k d i s s + k i n t ) M b   . (12)

The following differential equation represents the temporal variation of the molar concentration of remaining binding sites on the cell surface on which particles can still attach. It is proportional to dissociation and the internalization of bound particles and association rate of binding sites and unbound particles and also the rate of change of internalized particle related to internalization constant, respectively, as follows: ∂ M b s ∂ t = ( k d i s s + k i n t ) M b − k a s s M b s M ϵ ,   a n d   ∂ M i ∂ t = k i n t M b . (13)

The boundary conditions for the initial and at a later time t are as follows: initially, when all particles are present on the surface of the cell, the molar concentration of free, bound, and internalized particle bears no value, i.e., at t = 0, M (0,r) = M_b (0,r) = M_int (0,r) = 0 for 0   ≤ r < R ; and at a later time t, the molar concentration of particles free in spheroid volume holds the same value as the molar concentration of nanoparticles outside the spheroid   ( M 0 ) multiplied with porosity of spheroid, i.e., t = t, r = R, we have M(t, R) = M₀ ϵ (R), with no radial variation of unbound particle intracellularly   ∂ ∂ r ( M ϵ ) ( t , 0 ) = 0. Here, M₀ is considered to be uniform and equal to the averaged concentration of particles due to mixing in stirred vessels under the experimental conditions. Solving the above set of partial differential equations, we get the total number of particles that are retained in spheroid at time t as   M t ( t , r ) = p ( r ) [ M ( t , r ) + M b ( t , r ) + M i n t ( t , r ) ] ,   where p(r) is the piecewise linear window function. p ( r ) is used to consider the particles near the spheroid’s outer rim that need to be removed with the value of r 1     approximated to   0.95 R . It acquires unity when all the particles are internalized and less than unity as per the variation of their radial position from the spheroid surface and center, i.e., p ( r ) = 1 for   r < r 1 ; p ( r ) =   1 −   r − r 1 R − r 1   for   r 1 < r     ≤   R .

Once the nanoparticles are internalized, diffusion sets are limited to the targeted delivery. Diffusion plays an indispensable role in how the nanomedicine gets transported to multilayers from cellular scale to extracellular scale. Diffusion is driven by gradients in concentration. It involves the penetration of nanomedicine to multilayers of cells prior to making it to the center. Out-turns of diffusion barriers in multilayer tissues such as cellular compaction, efflux system, and gap junctions have been solved via various models (Costa et al., 2016). The model by Gao et al. (Gao et al., 2013), for the diffusion coefficient of intercellular, porous spaces into the spheroid is represented in three steps. Initially, free diffusion D 0 nanoparticles in water at 37 °C are calculated using the Stokes–Einstein equation. This coefficient at a specified temperature is given by D 0 = k B T 6 π η a , where k_B is the Boltzmann constant, T is the temperature in kelvin,   η is the viscosity of water, and a is the radius of the particle. Now, this is followed by diffusion in the interstitium   (   D i n t ) . It represents the diffusion in a porous gel matrix, which acts as an extensive connective network for cell-cell interactions and is responsible for the formation of an extracellular matrix via matrix protein. It depends on the ratio of particle radius ( a ) to pore radius ( a f ) and the square root of volume fraction of tumor interstitium matrix ( ϕ i n t ) . The corresponding relation can be expressed as D i n t = D 0 e ( −   a a f ϕ i n t ) , where the volume fraction of tumor interstitium matrix ( ϕ i n t ) is defined as the ratio of the volume occupied by the particle diffused in the ECM to the total volume of all particles present before diffusion and ϕ i n t > 0 always. This ratio is always less than equal to unity depending on the volume of particles that are diffused. For the calculation of   ϕ i n t , the value of collagen is used, which can be obtained by the multiplication of interstitial collagen concentration and its effective cell volume. For movements of nanoparticles to occur through the ECM, the typical mesh size of the tumor matrix should be comparable or larger than the size of the nanoparticles. Hence, the ratio   a a f   is usually less than unity. Collectively, this will lead to the overall value of ( a a f ϕ i n t ) less than unity. Now, incorporating diffusion at interstitial cellular scale, which mainly depends on the cellular density ϕ and porous spaces inbetween the cell ( 1 − ϕ ) can be represented as D i = D i n t ( 1 − ϕ ) 2 ; i.e., D i = D i n t   ϵ 2 , where D_i is the interstitial diffusivity constant for diffusion in a porous gel matrix with immobilized cells and ϵ is the porosity in spheroids and defined as the ratio of the volume of immobilized cells to the total volume of cells in the tumor matrix and lies in the range of 0 < ϵ < 1 . The interstitial diffusivity ( D i ) depends directly on the porosity factor, which acquires values less than unity, resulting in slightly lesser diffusion at the cellular scale. Thus, this entire process of diffusion can be summarized with an inequality relation: D 0 > D i n t > D i , evidently confirming the idea of the uneven rate of particle diffusion till they reach the spheroid’s core.

Taking the shape factor of spheroid into account consisting hindrances in its pores, steric hindrances due to the presence of ligand, and the tortuosity due to structural non-uniformity, the mathematical model for the nanomedicine diffusion has been given by Goodman et al. (Goodman et al., 2008). Similar to the model by Gao et al., initially, this model also takes free diffusion in the unbound medium that is not mediated by any force and potential due to cell-cell and cell-ECM interaction into consideration. The free diffusion given by the Stokes–Einstein relation can be written as   D 0 = k B T 6 π η a   , where D₀ is the diffusion coefficient in the unbounded medium, k_B is the Boltzmann constant, η is the viscosity of the liquid, and T is the absolute temperature. Now, depending on tortuosity τ ( ϵ ) , shape factor F, hindrances due to hydrodynamics, and the steric effects of the diffusion coefficient in the porous media L ( λ ) , the corresponding relation for effective diffusion constant can be presented as D   =   D 0   L ( λ ) F τ ( ϵ ) . Here, τ ( ϵ ) is related to the curvature consisting of structural aspects of the spheroid and has a profound impact on the diffusion and hydrodynamics in porous media. It is responsible for the increase in diffusion path length in the spheroid and related to the volumetric porosity ( ϵ ) as 1 τ ( ϵ ) = 1 − 2 3 ( 1 + ϵ ) ( 1 − ϵ ) 2 3 , where ϵ is the ratio of the volume of immobilized cells to the total volume of cells in a tumor and bears a value in the range 0 < ϵ < 1 . As defined, the inverse proportionality between them leads to τ ( ϵ ) >1. Mobility of particles in spheroid   λ is defined as the ratio of particle’s radius   a to the pore radius r p ; that is, λ = a r p . The typical size of pore radius is comparable or greater than the particle size as a < r p ; therefore, λ < 1 . The ideal situation of free diffusion is often affected by hindrances in the intercellular medium and tumor spheroid structure. Steric hindrances are defined as the obstructions in the path of a particle due to the presence of surrounding particles and ligands, often slowing down its motion and may stagnate the diffusion in the pore. This steric reduction along with the effects of hydrodynamic forces can be expressed as L ( λ ) = ( 1 − λ ) 2 ( 1 − 2.1004 λ   +   2.089 λ 3 − 0.948 λ 5 ) . The structure factor (F) incorporates the obstructions in the spheroid pore and holds a value F > 1. Thus, the combined effect of all these defines the demeanor of diffusion in spheroid, stating an overall lesser value of D than D 0   upholding the real phenomenon.

Binding sites are the receptors on the surface of cells that help the drug molecules and ligands to bind with them for better signal transduction pathways. One of the most crucial parts of targeted drug delivery is the appropriate availability of these binding sites on the surface of cells with which drug molecules can attach and eventually internalize. For impaired efficacy of drug delivery to the target, calculation of the concentration of available binding sites on the cell surface is a requisite. For the porous media, a parallel pore model was developed, which established a relationship among the molar concentrations of the available binding sites on the cell surface ( M b s ) taken initially at t = 0 with the spheroid structure, porosity, structural dimension of a particle, and pore radius as (Goodman et al., 2008)   M b s ( 0 , r ) = 2 π k β β ϵ a 2 r p N A , where a is the radius of the particle, r_p is the pore radius, N_A is the Avogadro’s number, and the density of remaining binding sites present on the cell surface is represented by β , whereas k_ꞵ is related to the variation in binding cell density on the monolayer cell culture surface to the cells in spheroid.

The surface of the tumor is well decorated with binding sites as receptors for the targeting ligand. Once the drug molecules bind to its surface, the next concern is about the calculation of the total number of cohesive particles that are bound and internalized. A model with several assumptions for a single cell considering bound and internalized particles to determine the number of cohesive particles, rate constants, and the number of binding sites is used. The concentration of particles surrounding the cell and the mean concentration of particles in the vessel are equal. A cell’s surface area is not inhibited by any other cells. Due to the continuous turnover of the cell membrane, the cell regenerates potential binding sites as the particle internalizes, as a result of which binding sites are taken to be constant, and any exocytosis will end up in the reduced value of k_int (Goodman et al., 2008). The model provides the relations representing the variation of number of bound particles with time proportional to the decreasing rate of dissociation and internalized particle and the positive influence of associated binding site on the cell surface along with total particle outside the spheroid as follows: d N b d t = − ( k d i s s + k i n t ) N b + k a s s M 0 N b s . (14)

The remaining number binding site N_bs on the cell surface is the difference between the number of available binding sites S, which were present initially, and the number of bound particles on the cell surface N_b: N b s =   S − N b   ,   a n d     d N i n t d t = k i n t N b   , (15) where N_int is the number of particles that are internalized and proportional to the internalization rate constant of bound particles. The initial number of binding sites (S) can be related to the effective segment of cell surface area that is available for binding ( β ) as β = S a 2 4 R c 2 , where   R c is the individual cell’s radius, and a is the particle radius. On the other hand, the growth of the number of bound particles and the total number of adhered particles can be given by solving Eq. 15. Thus, we get the number of bound particles ( N b ) and the total number of adhered particles ( N t   ) : N b = α ( 1 − e − γ t ) , N t ( t ) = N b ( t ) + N i n t ( t ) = α [ k i n t t   + γ − k i n t γ ( 1 − e − γ t ) ] , (16) where α and γ are defined as follows: α = k a s s M 0 S γ ,   a n d     γ = k a s s M 0   + k d i s s + k i n t   . (17)

Here, as γ is the sum of all the rate constants,   γ > k a s s and γ > k i n t individually. Also, S and M 0   are definite quantities and greater than unity. Hence, α < 1 , which leads to the fact that N b   a n d   N t   are less than unity and therefore sustaining our aim of greater internalization of particles into the spheroid. Solving Eq. 16 and using experimental data, the value of particle binding to the cell, rate constants, and the number of binding sites can be calculated. Thus, the mathematical models provide us the tool and allow us to calculate all the parameters such as concentration of nanoparticles, binding sites, diffusion coefficient at distinct scales, number of adhered particles, rate coefficients, and number of binding sites which are needed for successful targeted drug delivery. This quantitative approach to the problem makes us well equipped to deal with the system qualitatively as well by studying the physicochemical properties of drug carriers.

Quantitative analysis of nanomedicine in all the physical and mathematical regimes leads the way for the study of physicochemical features of the drug. The study of these properties of nanomedicine can help in establishing a suitable platform for target-oriented delivery of drugs with more adaptability (Liu et al., 2012). The following section will illustrate the key features of nanomedicines aiming at providing efficient drug delivery mechanisms, attaining the objectives, and enhancing bioactivity (Table 3).