Modular Representation of Physiologically Based Pharmacokinetic Models: Nanoparticle Delivery to Solid Tumors in Mice as an Example

: Here we describe a toolkit for presenting physiologically based pharmacokinetic (PBPK) models in a modular graphical view in the BioUML platform. Firstly, we demonstrate the BioUML capabilities for PBPK modeling tested on an existing model of nanoparticles delivery to solid tumors in mice. Secondly, we provide guidance on the conversion of the PBPK model code from a text modeling language like Berkeley Madonna to a visual modular diagram in the BioUML. We give step-by-step explanations of the model transformation and demonstrate that simulation results from the original model are exactly the same as numerical results obtained for the transformed model. The main advantage of the proposed approach is its clarity and ease of perception. Additionally, the modular representation serves as a simpliﬁed and convenient base for in silico investigation of the model and reduces the risk of technical errors during its reuse and extension by concomitant biochemical processes. In summary, this article demonstrates that BioUML can be used as an alternative and robust tool for PBPK modeling.


Introduction
The use of mathematical modeling methods has become an integral part of drug discovery and development [1]. Pharmacokinetic (PK) models can be applied to calculate the effective half-life of drugs after intravenous or oral administration to humans or animals, adjust the initial dose, and determine the alternative doses or dosing regimens in the target population.
Two main approaches can be used to analyze PK data [2]: • noncompartmental analysis estimates the exposure to a drug (such parameters as the area under the plasma drug concentration-time curve, maximum drug concentration, and time at maximum concentration) directly from observed individual profiles. • model-based compartmental analysis predicts the concentration-time curve using kinetic models.
The conceptual basis of the noncompartmental approach is that no assumptions are made about the underlying structure of the model. The main focus of the analysis is on the assessment of drug clearance and volume of distribution, rather than the mechanisms of drug disposition. The clearances and volumes calculated are numerically identical for both approaches. The major difference relates to the parameter used to describe how long a drug remains in the body. In the first approach, the mean residence time is calculated as the time that bisects the area under the plasma concentration-time curve. In the second approach, the elimination half-life (the time it takes for 50% of the administered drug to be eliminated from the body) is used [3].
The simplest form of the compartmental PK model is defined by one or two compartments [4]. The fundamental assumption behind the one-compartment model (representing the whole body) is that equilibrium is rapidly achieved within the body after the drug administration. In this case, the concentration profile is expressed as an exponential function of time. The two-compartment model (the first compartment is considered as a lumped sum of rapidly-perfused tissues, such as the liver, kidney, heart, and lung, while the second one represents a lumped sum of slowly-perfused tissues, such as the muscle, fat, and skin) can be used to describe two main phases, including the distribution and the elimination of drugs, which are also characterized by exponential dynamics.
More complex physiologically-based pharmacokinetic (PBPK) models consist of several compartments, each of which represents a specific part of the body: a group of organs or tissues (for example, tissue with slow perfusion), a single organ or tissue (e.g., liver, lungs, etc.) or some area of an organ or tissue (for example, intracellular space) [5,6]. These models are built using differential equations parameterized with known physiological variables, which can describe the absorption, distribution, metabolism, and excretion of drugs [1].
PBPK models are typically stipulated as a set of interconnected vessels with ideal mixing, where both biochemical and transport mechanisms are given as black-box, empirical relations. The kinetics of novel drugs can be studied more systematically with mechanistic biochemical models in entire organisms [7].
An application of molecular modeling techniques including molecular docking in the combination with PBPK models and high throughput experimental or clinical data is crucial in the process of drug discovery nowadays. The molecular docking approach, in particular, enables predicting the detailed interaction pattern between a chemical compound and its protein target or protein and protein, thereby providing the kinetic parameters for the PBPK model. The impact of the parameter values on the PBPK modeling results may be evaluated via sensitivity analysis. However, the molecular docking approach might not be accurate, and more advanced binding free energy calculation methods should be applied [8,9].
One of the applications for PBPK models is the study of nanoparticle delivery to tumors. Cancer is a major public health problem with millions of new cases each year resulting in hundreds of thousands of deaths. In the past three decades, significant efforts have been focused on the improvement of the targeted delivery of new highly effective therapeutic agents for the growing tumor. The development of multifunctional nanoparticles that can passively penetrate into tumor tissues without affecting healthy ones provides the most promising approach to solving the issue [10]. Despite these supposed benefits, the introduction of nanoparticles into clinical practice is still limited due to the low efficiency of delivering the administered dose to the tumor (0.70% according to the meta-analysis) [11].
Studying the biodistribution and interaction of engineered nanoparticles with cancer cells, as well as evaluating the effectiveness of treatment using these nanoparticles, covers a wide range of experimental approaches. Current technologies provide the following experimental techniques: • a combination of 2D and 3D optical tomography for in vivo imaging. • methods of fluorescent and traditional histology, spectroscopy, flow cytometry, and electron microscopy. • methods for obtaining transcriptomic and proteomic data to assess the expression of receptors for various compounds (sugars, lipids, albumin, etc.).
The application of PBPK models in cancer nanomedicine is only at a premature stage to date [17,20]. However, in recent years, a number of models have already been created to study the delivery of properties of various types of nanoparticles in tumor-bearing mice [11,21,22]. The advantages of such nanomedicine-oriented PBPK models are that [20]: • they allow studying the distribution of nanoformulations at the site of their action in the body. • they provide tools to better understand the microenvironment complexity, phenotypic diversity, and genetic heterogeneity of tumors. • they may prove useful in nanomedicine development to achieve optimum nanoformulation concentration profiles needed in tumor-targeted cells. • they could directly be exploited for specific cancer patient groups (e.g., pediatric population) to help solve issues of specific drug development and administration dosage scheme protocols.
In this study, we consider a generic PBPK model describing nanoparticle distribution in tumor-bearing mice [11]. This model characterizes endocytic uptake in the liver, lungs, spleen, kidneys, muscles, and other tissues of the mice. The biodistribution of nanoparticles in each compartment is also regulated by the ability to cross the capillary membrane of the organs [19]. As an additional compartment, this model contains cancer cells that are also capable to absorb and accumulate nanoparticles.
The authors of the study provided a model as the code for the Berkeley Madonna platform (https://berkeley-madonna.myshopify.com/ accessed on 1 April 2022), which is not intended for visualization of biochemical systems. Such a representation complicates the understanding and analysis of the model for experts in biomedicine, and also complicates its use and development for mathematical modelers.
Herein, we describe the consistent transformation of the model [11] to a modular view using the Biological Universal Modelling Language (BioUML) platform, an integrated web framework for systems biology that provides a wide range of capabilities, including access to biological databases and tools for visual modeling, simulation, parameters fitting and model analysis [23]. The procedure comprises the transition from a system of differentialalgebraic equations to a chain of biochemical reactions and translation of the common compartmental structure into the form of separate modules. The visually modified version of the model reproduces the simulation results presented in the original study confirming the correct transformation of the model performed in BioUML.
We also point out that the visual modeling of PBPK is a novelty in the research field. Previously, text-based software packages such as acslX, Berkeley Madonna, MATLAB, and R language were mostly used to build this type of model [24]. While tools with visual modeling support simplify the interpretation, edition, and extension of PBPK models reducing the risk of technical errors. In addition, PBPK models reconstructed in BioUML are consistent with existing modeling standards and can be exported to the COMBINE archive [25].
The obtained results are methodological (platform) prerequisites to effectively develop new technologies for targeted and personalized delivery of nanotherapeutic agents in the treatment of cancer (for example, breast cancer and prostate cancer). The proposed modular model will be used to study the mechanisms of the substance transport in tumors and the delivery of drugs to them, taking into account the processes of nanoparticle biodistribution depending on their material, as well as the growth and metabolism of the tumor.

Structure of the Nanoparticle Delivery Model
The model of nanoparticle delivery to solid tumors ( Figure 1) [11] includes 10 compartments: arterial and venous plasma, lungs, spleen, liver, kidneys, brain, muscle, tumor, and remaining tissues. Besides plasma and brain, all compartments are divided into three subcompartments: capillary blood, tissue interstitium, and endocytic/phagocytic cells or tumor cells. The diagram in Figure 1 demonstrates that after intravenous administration, nanoparticles enter the capillary blood of the lungs with the bloodstream then move to the arterial plasma and disperse to all other organs of the mice with excretion via both biliary and urinary pathways. The model takes into account the uneven distribution of injected nanoparticles between capillary blood and organ tissue, membrane-limited transcapillary transport, as well as nonlinear endocytic uptake, and first-order exocytic release.
The rate of intravenous administration of nanoparticles is calculated by the injection time (denoted by Timeiv), the injected dose (PDOSEiv), and the bodyweight (BW) according to the formula (all parameters in the current section have the same notation as in the original model, and are listed in Table 1 for convenience, including values and units): otherwise. Volume of the capillary blood in the lungs 0.00007 L VLuC Lung volume (fraction of the bodyweight) 0.007 -VLut Volume of the lung interstitium 0.00007 L Further distribution of nanoparticles over compartments is determined by the rate of blood flow and metabolic processes within them. Next, we consider the modeling of these processes using the lungs as an example.
One of the main parameters of the lung compartment is its volume VLu, calculated as a fraction VLuC = 0.007 of the weight: It is assumed that the value of VLu is equally divided by the volume of capillary blood VLub and the volume of tissue interstitium VLut: Denoting the number of nanoparticles in these subcompartments by the parameters ALub and ALut, respectively, we obtain formulas for calculating the concentrations:

Now let
ALuRES be the amount of nanoparticles being taken up from the lung interstitium by phagocytosis. Then the total amount of nanoparticles in the lung compartment ALung is defined by the sum: Describe the reactions that affect the values ALub, ALut, and ALuRES. If CV is the concentration of nanoparticles in venous plasma and QC is cardiac output, then the rate of nanoparticle entry into the capillary blood of the lungs is determined by the product of these quantities (reaction r1 in Figure 2): The next step in the distribution of nanoparticles is diffusion through the capillary wall into the tissue interstitium. This process depends on the membrane-limited permeability coefficient, PALuC (r2in, Figure 2 while the rate of the reverse process is calculated by the formula (r2out, Figure 2): where PLu is the tissue:plasma distribution coefficient of nanoparticles. The endocytic/phagocytic uptake is described by the Hill function: with parameters KLuRESmax, KLuRES50, and KLuRESn, whereas the release from endocytic/phagocytic cells is determined by a constant coefficient KLuRESrelease. The rates of the corresponding reactions are given by the equations (r3in and r3out, Figure 2 The last reaction to consider is the outflow of nanoparticles from capillary blood into arterial plasma, which occurs at the rate (r4, Figure 2): Taking into account the rate equations above, we derive the system of ordinary differential equations representing nanoparticles distribution between subcompartments of the lungs: The implementation of these formulas in Berkeley Madonna, provided by the authors in the article [11], is shown in Figure 3. Despite the fact that this implementation contains comments describing parameters and expressions, it is difficult to understand the sense of the entire model. However, as follows from representation (6), the right side of the major equations provides a way to consider it as the sum of individual reaction rates. The representation makes the model view more explicit using the capabilities of the visual modeling in BioUML [23].
The next section includes a brief description of the methods that we used to convert the original model [11] to a modular form with a visual representation of the system of equations as a chain of biochemical reactions.

Visual Modeling
Representation of the biological system as a graphical diagram using a tool supporting visual modeling can significantly simplify the procedure of the mathematical model construction. According to the visual modeling paradigm, each component of the diagram corresponds to a certain mathematical object (variable, reaction, equation, etc.). A user can edit the properties of each element on the diagram: for instance, set initial values of the model variables or specify kinetic laws of enzymatic reactions. The BioUML platform (https://sirius-web.org/ accessed on 1 April 2022) [23] is based on the visual representation and certain properties of the diagram elements automatically generate a highly optimized Java code that is employed for numerical simulation of the model. Depending on the model type, BioUML provides a list of available simulation engines. For example, reaction pathways can be simulated as a system of ordinary differential equations (ODEs) or as a Gillespie-type stochastic model. Each simulation engine provides a list of suitable solvers, in particular, solvers for systems of ODEs include JVODE (a package CVODE ported from C to Java [26], which uses multistep Adams-Moulton method), RADAU5 solver [27], as well as classic algorithms (Euler, Dormand-Prince [28]).
An essential component of visual modeling is a graphical notation that enables formalizing and fully building the model of a biological process or phenomenon. The visual model can be reconstructed using several types of diagrams providing opportunities to describe different aspects of both structure and function of the complex biological systems considering diverse levels of their details.
As part of the visual modeling development in BioUML, an extension of the Systems Biology Graphical Notation (SBGN) standard [29] was suggested for comprehensive graphical representation of Systems Biology Markup Language (SBML) models [30] taking into account equations, discrete events, math functions, and constraints [31] (Table 2). There is also implementation to integrate table data into the model which may determine dynamics changes of variables. For example, a report of the treatment introduction and drug dose can be set as a table for pharmacokinetic models.

Equation
Mathematical equations of the model: Discrete event Events describe an instant change of the model variable when a specified condition is satisfied, e.g., the change of medication rate after a certain time point.

Function
Function receives argument values, calculates, and returns the result of the given expression.

Constraint
A condition set in the constraint is checked during the simulation. Depending on solver settings a simulation may be stopped or an error message will be shown if the condition is violated.

Interface ports
There are three types of interface ports in the BioUML to support the visual representation of modular models (see subsection "Modular modeling" below): contact (gray), input (green), and output (red) ports.

Modular Modeling
Modularity, as a principle of the complex biological systems organization [32], has served as a basis for the creation and intensive development of the modular approach to simulate such systems over the last two decades [33,34]. In the framework of this approach, a system is considered a set of interconnected subsystems (modules). Each module is a mathematical model and diverse numerical methods (simulation, parameter fitting, sensitivity analysis, etc.) can be independently applied to its analysis. Integration of these modules enables building a more complex modular model of the whole system. Modules may have different mathematical formalisms (discrete, deterministic, stochastic systems) in order to describe internal processes with different details' extent of their mechanisms.
A special type of modular diagram has been suggested for their visual modeling in BioUML. This type enables specifying connections between modules. The subset of variables used for that is called the module interface and is represented by the corresponding ports [31] (Table 3): • input port-mathematical variable used in the module, but calculated outside; • output port-mathematical variable calculated in the module and may be used outside; • contact port-mathematical variable that can be modified both inside and outside the module.
Besides the ports listed above, modular diagrams also include all math elements and tabular data represented in Table 2. Table 3. Designation of the model elements on modular diagrams in the BioUML. Adapted from https://doi.org/10.3390/ijms221910353 accessed on 1 April 2022.

Module
Module encapsulates the mathematical model (submodel) of a particular subsystem as an element of the hierarchical structure of the more complex model of the whole system.

Interface ports
Ports represent the interface of the module through which it can be connected with other modules or with the hierarchical model itself.
The color of the port defines its type (the same as in Table 2).

Connection
Connections can be established between ports to aggregate modules with each other. Directed connections are established between output and input ports, while undirected connections are between contact ports.

Bus
Buses are auxiliary elements of the modular diagram that can be used as intermediate points in connections. They allow decreasing the number of intersections between connections and make the structure of the model clearer.

Antimony-Textual Representation of the Model
Antimony is a tool for numerical model construction as a convenient human-readable text format that supports most of the SBML features ( Figure 4) [35].  SBGN. (B). The corresponding Antimony code which includes the definition of compartments (Capillary_blood and Tissue_interstitium), species (amounts of nanomaterial in capillary blood, Alub, and in tissue interstitium, Alut), parameters (PALuC and QC), reactions (R1) with the related kinetic laws (a math equation to the right of the reaction declaration), and properties for the visual model representation, including titles of the model elements, and specification of the SBGN notation.
An extension of the Antimony for specification of reaction components as well as for some other properties of SBGN diagrams has been implemented in the BioUML according to the next rule: @entity_id.property = value where entity_id-an identification of the particulate model element, property-property of the element, and value-value of the property. The proposed extension is similar to the annotation principle in SBML format where SBML supporting software can store any auxiliary information in an XML file.
Currently, the BioUML supports the next structural properties for the annotation in text format: • sbgnType-defines the SBGN entity type (macromolecule, nucleic acid, simplex chemical or complex, unspecified); • sbgnViewTitle-defines additional properties of an SBGN entity (e.g., multimer). If an entity is complex, it also defines the inner elements. The Transpath database specification is used to denote entities and complexes in the BioUML (https://genexplain.com/ wp-content/uploads/2017/04/TRANSPATH-Documentation_2012.2.pdf accessed on 1 April 2022). For instance: A:B-complex comprising two entities, A and B; A{p}-entity A in the phosphorylated state; (A)3-multimer A consisting of three subunits.
A precise representation of the model for a biological system as a graphical diagram considering extended functionality for SBGN notation or as a code using Antimony depends on the context and aim of the modeling. The BioUML platform supports both types of representation which are built in the platform and unambiguously correspond to each other meaning that a change in Antimony code leads to automatic modification of the diagram and vice versa.

Reconstruction of the Computational Modular Model
According to the modular approach [36], a biochemical system is considered a set of functioning interconnected subsystems. Such visualization adequately reflects the structure of biological objects at different levels of their organization from cell to organism [37,38]. PBPK models are usually parameterized taking into account the known physiology of the studied organisms [1]. The structure of the models is divided into compartments (organs and tissues in the body) connected by the circulating blood system and, therefore, can be easily represented using modular mathematical formalism. In this case, the modeling of individual compartments is based on the same differential-algebraic equations, with different values of the parameters. Therefore, a logical step is isolating the common equations as a separate module that receives input values specific to each compartment. To avoid errors in the reconstruction of the model [11] and obtain identical simulation results, we performed the model transformation in several steps: • transition to a modular structure. • transition from a system of differential-algebraic equations to a chain of biochemical reactions. • finding the common structure of compartments as a separate module Below we describe the transformation procedure step-by-step in more detail.

Transition to a Modular Structure
To import the model [11] into BioUML, we used the Antimony language [35], supported by the platform and intended for textual definition of biochemical systems. The Antimony and Berkeley Madonna functions have a similar syntax, which, for example, is demonstrated by the formulas for intravenous administration of nanoparticles in Figure 5. Translation of the model code into the Antimony language, divided into several parts corresponding to individual compartments, gave us a set of computational submodels, each of which was used as a separate module. Since BioUML supports Antimony integration with a diagram view of models, these submodels were automatically converted into visual diagrams ( Figure 6). The intermediate version of the modular model, which we obtained by connecting the common variables of the modules, is similar to the final version ( Figure 7B), except that the nanoparticle dosing formulas ( Figure 5), initially placed into a separate module, were subsequently included into the module "Venous plasma" to fully comply with the scheme in Figure 1. Since this intermediate version is not essential, we omit its visualization in the current article and suggest referring to the web implementation of the model for details: https://sirius-web.org/bioumlweb/#de=data/Collaboration%20(git) /Nanoparticles/Diagrams/02/PBPK%20model accessed on 1 April 2022.

Transition from a System of Differential-Algebraic Equations to a Chain of Biochemical Reactions
As mentioned above, the change in the amount of nanoparticles in different compartments and subcompartments of the model is based on the summation of local rates of their distribution. The system of Equation (6) with the amount of nanoparticles in venous (AV) and arterial (AA) plasma can be defined as a chain of reactions with the rates v 1 , v in 2 , v out 2 , v in 3 , v out 3 , and v 4 (1)- (5): In such formulation, biochemical processes can be specified in the SBML format [30], as well as visualized using the graphical notation SBGN [29] (Figure 8). Figure 8. The lung compartment in BioUML using Systems Biology Markup Language (SBML) and SBGN standards and Antimony code. Note that the "sbgnType" tag is not included in the Antimony syntax, as well as specifying the different types of ports (>BW, <QC). We added these extensions to the Antimony implementation to support SBGN objects and allow for modular modeling.
Note that due to the change in the module's structure, the general scheme of the entire model was also changed. While in the first version of the model, the concentrations of nanoparticles in arterial and venous plasma were transferred through directed connections, in the second version, these connections are undirected and correspond to the amount of nanoparticles in the plasma with the parallel transfer of plasma volumes to calculate concentrations directly in the compartments.

Finding the Common Structure of Compartments as a Separate Module
The formulas of membrane-limited transcapillary transport, endocytic uptake, and exocytic release of nanoparticles in different compartments have the same form (2)-(4), but different values of the kinetic constants. The distribution of nanoparticles between capillary blood, tissue interstitium, and endocytic/phagocytic (or tumor) cells for most compartments occurs according to the same scheme as for the lungs (Figure 8), with a difference in the incoming (arterial instead of venous) blood flow. Therefore, extracting the reactions of nanoparticle transport between subcompartments into a separate module ( Figure 9A) and defining input ports for specific parameters and variables in this module results in a common structure module that can be used to model individual compartments. Using this structure, we designed the lungs, spleen, kidneys, brain, muscle, tumor, and the compartment responsible for the rest of the body ( Figure 9B-H). This step resulted in the following differences: • inflow and outflow for the lungs are determined by the amount of nanoparticles in the venous and arterial plasma, respectively, while for the remaining compartments, these flows are determined vice versa; • to convert the amount of nanoparticles in plasma into their concentration and calculate rates of the form (1), we pass the value of plasma volume PV to the compartments of the model: we consider venous plasma volume VPv ( Figure 7B) for the lungs and arterial plasma volume APv for other compartments; • input ports on the left and right side of the common structure module shown in Figure 9B-H correspond to the individual compartment constants (with the same notations as in Table 1 for the lungs); since in the model [11] the brain compartment does not contain the subcompartment of endocytic/phagocytic cells, the parameters associated with it are not transferred and their default values remain zero • in addition to the common equations, the kidney compartment includes the reaction of nanoparticles excretion through the urinary tract ( Figure 9D), proceeding at a rate: where AKb, VKb, and BW denote parameters, the values of which come from the common structure module and determine the amount of nanoparticles in the capillary blood of the kidneys, the volume of capillary blood, and bodyweight, respectively. KurineC is the rate constant of nanoparticles excretion through the urinary tract.
• the tumor compartment contains formulas for calculating the amount (ATumort) and concentration (CTumort) of nanoparticles in the tissue as a percentage of the administered dose (DOSEiv): CTumort DOSEiv where 0.001 is the coefficient for converting tissue volume units from liters to grams.
• to use the cardiac output (equal to the blood flow in the lungs) and bodyweight outside the compartments shown in Figure 9, we provided for the transfer of these values from the common structure module, where they are determined, to the top level of the model (see modular representation in Figure 7B) via the lung compartment ( Figure 9B, port Q at the bottom of the compartment) and the rest of the body compartment ( Figure 9G, port BW at the top of the compartment), respectively.
Only the liver cannot be expressed through the common system of reactions. The difference with other compartments is that this organ contains Kupffer cells, which are able to phagocytize nanoparticles directly from the capillary blood [11]. The sequence of reactions, in this case, is shown in Figure 10A. In addition to the processes r1, r2in, r2out, r3in, r3out, and r4 with equations similar to (1)- (5), this sequence includes the excretion of nanoparticles with bile (r5): where ALt is the amount of nanoparticles in the liver tissue, VLt is the tissue volume, BW is the body weight, and KbileC is the rate constant of nanoparticles excretion with bile.  Table 1, corrected for the corresponding module.
To complete the description of the model, we should note the modules of venous and arterial plasma, which are presented in Figure 10B,C, respectively. The final version of the modular computational model is available at: https://sirius-web.org/bioumlweb/#de= data/Collaboration%20(git)/Nanoparticles/Diagrams/04/PBPK%20model accessed on 1 April 2022.

Testing the Model after the Modification
We benchmarked the modular version of the model using the next test cases: the case of the healthy mice in which the tumor compartment is omitted, and the case of tumorbearing mice with the parameters of the tumor compartment given in Table 4. The model dynamics in the first case were calibrated by the authors using the biodistribution data for venous plasma, lungs, liver, spleen, and kidneys in healthy mice intravenously injected with 13 nm gold (Au) nanoparticles (AuNPs) [39]. In the second case, the parameters of the tumor compartment were originally identified using experimental data on tumor suppression in mice by photoacoustic treatment with gold nanorods (AuNRs) of 9-11 nm in diameter [40]. Table 4. Values of the model parameters to simulate cases of distribution of the gold nanoparticles (AuNPs, 13 nm) in healthy mice and of the gold nanorods (AuNRs, 9-11 nm) in tumor-bearing mice.  Figure 11 demonstrates an exact quantitative agreement between the simulation results of the original model [11] and the final modular model in both cases. As was observed in the original study by Cheng and coauthors, the model simulations show a gradual decrease in the concentrations of AuNPs in the venous plasma, lungs, and kidneys as they accumulate in the liver and spleen of healthy mice over 168 h after administration. The dynamics of AuNRs in tumor-bearing mice are similar, with the difference that the concentrations of AuNRs in all organs are higher (due to the higher injected dose: 20 versus 0.85 mg/kg), while accumulation of AuNRs also occurs in the tumor tissue, thereby successfully reproducing the experimental dataset [40] for AuNRs biodistribution for up to 7 days after the administration. This in silico result confirms that nanoparticles can effectively target cancer cells and increase their accumulation in tumors. An interactive reproduction of these results in a Jupyter document (https://jupyter.org/ accessed on 1 April 2022.) is available at: https://gitlab.sirius-web.org/nanoparticles/models accessed on 1 April 2022. Figure 11. Simulation results of temporal changes in the concentration of the gold nanoparticles (AuNPs, 13 nm) in healthy mice and of the gold nanorods (AuNRs, 9-11 nm) in tumor-bearing mice, obtained for venous plasma, as well as lung, liver, spleen, kidney, and tumor tissues.

Conclusions
Herein, we describe in detail the sequential transformation of a physiologically based pharmacokinetic model of nanoparticles in mice with solid tumors to a modular computational form. The scheme of the original model included 10 compartments: arterial and venous plasma, lungs, spleen, liver, kidneys, brain, muscle, tumor, and remaining tissues. Therefore, the first logical step in the transformation was the definition of a separate computational module for each compartment and the indication of connections between modules in accordance with the movement of blood between organs and tissues of the body. All compartments in the original model, except plasma and brain, were initially divided into three subcompartments: capillary blood, tissue interstitium, and endocytic/phagocytic (or tumor) cells. The second step of the transformation was the transition from an explicit system of differential-algebraic equations to a chain of biochemical reactions and visualization of the resulting version of the model using the capabilities of the BioUML software. Finally, the nanoparticle transport equations in the original model had the same form for different organs, but a unique set of kinetic constants. Thus, the third step of the transformation was the allocation of the common structure of compartments into a separate module. Model versions obtained at each step reproduce the original simulation results and are available as SBML or COMBINE archive files in the web edition of BioUML (https://sirius-web.org/ accessed on 1 April 2022).
In summary, the present article focuses on the visual implementation of computational PBPK models using the BioUML software. We demonstrate that BioUML allows us to recode the published text-based model on the fly and simplify its representation for non-expert users based on the supported visual modeling and modular approach. This transformation of PBPK models actually enhances the ability of their application [24]. In this regard, we propose some thoughts to consider:

•
To share biological models (in particular, PBPK), it is important to make them available. The use of well-known biological standards is essential in this context. A translation of PBPK models from text-based formats (like acslX, Berkeley Madonna, MATLAB, and R language) to SBML and SBGN standards, which BioUML allows, makes these models accessible to a wide range of scientists and provides a quick start to work with them. Therefore, the presented description can be used as a guideline for the procedure. • Visual modeling serves as a convenient tool for exploring in silico models. This technique reduces the risk of technical errors during the model creation, use, and further development. In addition, the computational model in BioUML has the same view as the technical diagram of the model, which the developer can draw on paper to visualize the PBPK processes. This possibility is very useful for both beginners to learn PBPK modeling techniques and advanced developers. We believe that visual modeling is the future.

•
In this article, we not only demonstrate the BioUML capabilities for the PBPK modeling, but also transform the model of nanoparticle delivery to solid tumors in mice. In this respect, we consider the next possible roadmap for future research on this model: a. model simulations and analysis depending on the initial concentration of the administered nanoparticles dose. b. study the mechanisms of therapeutic nanoparticles transport taking into account their material and biodistribution depending on the tumor growth and metabolism.
This study was conducted as part of a project to develop new technologies for targeted and personalized delivery of nanotherapeutic agents in the treatment of cancer (https: //rscf.ru/en/project/21-75-30020/ accessed on 1 April 2022).
Funding: This research was funded by the Russian Science Foundation (project 21-75-30020).

Informed Consent Statement: Not applicable.
Data Availability Statement: The data are available on Gitlab at https://gitlab.sirius-web.org/ nanoparticles/models/ accessed on 1 April 2022.

Conflicts of Interest:
The authors declare no conflict of interest.