Mathematical and Computer Modeling as a Novel Approach for the Accelerated Development of New Inhalation and Intranasal Drug Delivery Systems

: This paper presents modern methods of mathematical modeling, which are widely used in the development of new inhalation and intranasal drugs, including those necessary for the treatment of socially signiﬁcant diseases, which include: tuberculosis, bronchial asthma, and mental and behavioral disorders. Based on the conducted studies, it was revealed that the methods of mathematical modeling used in the development of drugs are fragmented


Introduction
Socially significant diseases are those whose occurrence and spread are, to some extent, dependent on the socio-economic conditions of the population's lives. They also cause significant socio-economic damage due to a loss of productivity, treatment costs, disability, and population mortality.
According to the World Health Organization, noncommunicable, socially significant diseases (NCDs) include cardiovascular diseases (heart attacks and stroke), which are the cause of the majority of deaths from NCDs, and account for 17.9 million people per year, followed by cancer, with 9.3 million deaths per year, and chronic respiratory diseases (chronic obstructive pulmonary disease and asthma), with 4.1 million deaths per year [1]. The NCDs also include major mental and behavioral disorders. In addition to NCDs, infectious lung diseases such as tuberculosis, bronchitis and multidrug-resistant tuberculosis are included in the group of socially significant diseases. According to the statistics for 2021, about 10.6 million people fell ill with tuberculosis, of which 1.6 million people died [2]. Worldwide, tuberculosis is the thirteenth leading cause of death and the second leading cause of death after COVID-19 [2,3]; however, tuberculosis is curable and preventable. At the same time, multidrug-resistant tuberculosis poses a much more serious public health threat. In 2021, one in three people with drug-resistant tuberculosis sought medical care [2].
Thus, the development and introduction to the market of promising drugs based on modern delivery systems and new forms for the treatment of socially significant diseases is an urgent task today. The most promising technologies for drug delivery in the treatment of socially significant diseases are intranasal and inhalation technologies. Intranasal and inhalation delivery methods offer several advantages, including direct access of a drug to the systemic circulation, the potential for direct drug delivery to the brain through the olfactory nerves, convenience, and a rapid attainment of therapeutic drug levels in the drug to the systemic circulation, the potential for direct drug delivery to the brain through the olfactory nerves, convenience, and a rapid attainment of therapeutic drug levels in the bloodstream [4]. In addition, the use of nasal and inhaled forms of drugs avoids first-pass metabolism, allowing a drug to be delivered directly to the nasal cavity or lungs, thereby reducing the likelihood of adverse reactions from the gastrointestinal tract, and it removes the influence of pH and the fact of food intake on the rate of drug entry into the bloodstream [5]. The development of new pharmaceutical drugs, however, remains relatively slow due to several factors, including a low bioavailability, limited therapeutic efficacy, and a range of adverse reactions in the body, such as resistance, that may occur with high doses or repeated drug administration, as observed during clinical trials.
Mathematical modeling is recognized as an effective approach to expedite the development and improve the bioavailability of new pharmaceutical drugs. This approach involves the application of computational fluid dynamics methods, which enable the accurate prediction of the trajectory and deposition place of drug particles within the nasal cavity or lungs [6]. Furthermore, the development of cellular automaton models provides valuable insights into the release kinetics of drug substances from dosage forms [7]. Additionally, there are existing models [8] that comprehensively describe the pharmacokinetics and pharmacodynamics of these pharmaceutical drugs.
However, in the currently available scientific articles and literature reviews, there is no approach that combines the aforementioned methods of mathematical modeling. According to the authors, the development of a multiscale model that encompasses a unified set of mathematical modeling methods and incorporates sequential stages for predicting the trajectory and deposition sites of drug particles, as well as calculating the release kinetics of drugs from the dosage form and determining a drug's distribution throughout the body via blood flow, is a pressing task in the development of new drugs for the treatment of socially significant diseases.
The utilization of multiscale modeling as a novel approach in drug development offers a significant cost reduction in experimental studies. It enables the creation of a digital replica of a drug and facilitates the modeling of its structure and dissolution process in various environments, including the nasal cavity or lungs. Additionally, modern visualization methods for medical data, such as computed tomography (CT)- Figure 1 [9,10], allow for the precise digital representation of anatomical models of the nasal cavity or lungs. This, in turn, enables the accurate prediction of the flight path and deposition sites of drug particles during the developmental stage. Consequently, it provides the opportunity to optimize the composition of a medication to achieve the desired therapeutic effect, even prior to the commencement of clinical trials. The literature review conducted by the authors discusses the existing mathematical modeling methods extensively employed in the realm of drug development. The integration of these methods into a unified framework enables the development of a multiscale model capable of accelerating the process and market entry of novel drugs aimed at addressing socially significant diseases. The literature review conducted by the authors discusses the existing mathematical modeling methods extensively employed in the realm of drug development. The integration of these methods into a unified framework enables the development of a multiscale model capable of accelerating the process and market entry of novel drugs aimed at addressing socially significant diseases.

Mathematical Modeling of the Movement and Deposition of Drops/Particles during Inhalation and Intranasal Administration of a Medicinal Substance
Inhalation and intranasal routes of drug administration hold promise for the treatment of various diseases; however, achieving the desired therapeutic effect relies on the deposition of the drug substance in specific areas of the respiratory system [6,11]. Therefore, studying and predicting the deposition zones of medicinal substances administered via the intranasal or inhalation route is a crucial scientific and technical task to maximize the effectiveness of disease therapy.
Currently, in vivo studies provide reliable and accurate results regarding the movement and deposition of drugs in the respiratory system through different routes of administration. Modern methods of medical data visualization, such as scintigraphy or single photon emission computed tomography [9], are employed to investigate drug deposition zones. Radio-labeled particles enable the tracking of drug movement and determination of precise deposition sites; however, this approach is constrained by the high cost of imaging tools, complexities in particle labeling, and the potential risks associated with radioactive markers. Furthermore, in vivo studies are limited by patient-specific anatomical variations and the complexities of drug administration, making it challenging to analyze the acquired data and comprehend the mechanisms of drug transport and deposition in the respiratory tract [12].
An alternative to in vivo studies is in vitro diagnostics, which allows the study of drug deposition zones using simplified devices or physical replicas of the respiratory system [13]. This approach is cheaper and is not regulated by medical ethics; however, the in vitro approach is limited by the complexity of recreating an exact copy of the human respiratory system, taking into account the anatomical features [14].
With the increasing computational power in studying drug behavior in the respiratory system, in silico studies become viable, overcoming the limitations of in vivo and in vitro investigations. Mathematical models based on computational fluid dynamics equations have been employed to design inhalation and intranasal delivery devices [15,16] and to predict the movement and deposition of drug particles or droplets for various administration routes [17,18].
Computational fluid dynamics presents a promising approach in developing innovative inhalation and intranasal drug delivery systems; however, when utilizing this approach to study drug behavior, it is crucial to consider the numerous factors that influence the movement and deposition of particles or droplets in the respiratory system.

The Structure of the Respiratory System
The human respiratory system is a complex anatomical structure. It consists of the respiratory tract (i.e., the extrathoracic and tracheobronchial region) and lungs (i.e., the alveolar region) [19]. The airways, in turn, are divided into the upper and lower airways- Figure 2.
The respiratory system is divided into 23 segments (orders) from the trachea to the alveoli as shown in Figure 2. The upper respiratory tract includes the nasal cavity, the oral and nasal parts of the pharynx, and the larynx; the lower tract includes the trachea and bronchi (bronchial tree). Separately, it is worth highlighting the epiglottis, the function of which is to close the lumen of the larynx when swallowing.
The lower respiratory tract (tracheobronchial region) is the region of the respiratory tract from the trachea to the terminal bronchioles (16th order bronchi). The main function of the respiratory tract is to conduct air into the lungs, and to warm it and cleanse it of foreign particles and pathogens. The trachea is divided into two main bronchi, which, in turn, are divided into lobar (with two in the left and three in the right lung), and they are segmental up to the terminal bronchioles (i.e., a dichotomous type of branching). A feature of the structure of the trachea and large bronchi is the presence in their wall of cartilaginous semirings and rings, respectively. According to some studies [20], this feature leads to a greater deposition of particles in these parts of the respiratory tract. The respiratory system is divided into 23 segments (orders) from the trachea to the alveoli as shown in Figure 2. The upper respiratory tract includes the nasal cavity, the oral and nasal parts of the pharynx, and the larynx; the lower tract includes the trachea and bronchi (bronchial tree). Separately, it is worth highlighting the epiglottis, the function of which is to close the lumen of the larynx when swallowing.
The lower respiratory tract (tracheobronchial region) is the region of the respiratory tract from the trachea to the terminal bronchioles (16th order bronchi). The main function of the respiratory tract is to conduct air into the lungs, and to warm it and cleanse it of foreign particles and pathogens. The trachea is divided into two main bronchi, which, in turn, are divided into lobar (with two in the left and three in the right lung), and they are segmental up to the terminal bronchioles (i.e., a dichotomous type of branching). A feature of the structure of the trachea and large bronchi is the presence in their wall of cartilaginous semirings and rings, respectively. According to some studies [20], this feature leads to a greater deposition of particles in these parts of the respiratory tract.
The alveolar region (lungs) begins with the so-called respiratory bronchioles (bronchi of the 17th order) and connects through them with the final sections of the lower respiratory tract (terminal bronchioles). The respiratory bronchioles pass into the alveolar ducts, into which the alveoli open. The alveolar region consists of bronchi from 17 to 23 orders ( Figure 2). The main function of the alveolar region is gas exchange.
Modern methods of medical data visualization make it possible to build anatomically-similar models of the respiratory system and use them for in silico studies to predict the movement and deposition of particles or drops of drugs using computational fluid dynamics methods; however, at the moment there is a limited number of publications on the study of the movement of drugs in all areas of the respiratory system simultaneously due to the complexity of describing the nature of the air flow. It is noted in [10,21,22] that the reconstruction of a complete 3D model of the entire lung is currently not feasible due to the difficulty of segmenting the small airways using CT images; however, if the entire lung were fully segmented, it would be computationally-impossible to accurately model a complete tree of lungs, since this would require billions of mesh elements. Therefore, modeling of the particle flow through the lungs is carried out separately for the upper and for the lower respiratory tracts, as shown in Figure 3 (Stages 1 and 2). Modern methods of medical data visualization make it possible to build anatomicallysimilar models of the respiratory system and use them for in silico studies to predict the movement and deposition of particles or drops of drugs using computational fluid dynamics methods; however, at the moment there is a limited number of publications on the study of the movement of drugs in all areas of the respiratory system simultaneously due to the complexity of describing the nature of the air flow. It is noted in [10,21,22] that the reconstruction of a complete 3D model of the entire lung is currently not feasible due to the difficulty of segmenting the small airways using CT images; however, if the entire lung were fully segmented, it would be computationally-impossible to accurately model a complete tree of lungs, since this would require billions of mesh elements. Therefore, modeling of the particle flow through the lungs is carried out separately for the upper and for the lower respiratory tracts, as shown in Figure 3 (Stages 1 and 2).
The pulmonary trunk, that was presented in [10,21], ending before the beginning of the bronchi of the 16th order, had from 2 10 6 to 4 10 6 cells of the computational grid ( Figure 3a). However, in the developed complex acinar model in [22], spherical alveoli were attached to triple bifurcation units to represent generations of alveolar airways from 16 to 18 and 19 to 21 orders (Figure 3c,d), where the alveoli were separated by a minimum distance to allow for expansion by presenting septa that separate the alveoli. The final computational grid for single-generation, single-cell and multi-cell alveolar models contained 100,000 and 280,000 grid elements, respectively, while there were 2×10 6 grid elements in triple bifurcation units. Thus, the fragmentation of the lungs into the upper and lower airways, when calculating the flight and deposition of drug particles, is beneficial in terms of saving computational power and time. The pulmonary trunk, that was presented in [10,21], ending before the beginning of the bronchi of the 16th order, had from 2 10 6 to 4 10 6 cells of the computational grid ( Figure  3a). However, in the developed complex acinar model in [22], spherical alveoli were attached to triple bifurcation units to represent generations of alveolar airways from 16 to 18 and 19 to 21 orders (Figure 3c,d), where the alveoli were separated by a minimum distance to allow for expansion by presenting septa that separate the alveoli. The final computational grid for single-generation, single-cell and multi-cell alveolar models contained 100,000 and 280,000 grid elements, respectively, while there were 2×10 6 grid elements in triple bifurcation units. Thus, the fragmentation of the lungs into the upper and lower airways, when calculating the flight and deposition of drug particles, is beneficial in terms of saving computational power and time.

Computational Fluid Dynamics Method for Describing the Movement of Air in the Respiratory System
The use of computational fluid dynamics makes it possible to predict the properties of a liquid or gas flow by solving the differential equations of the basic conservation laws, presented in the form of the Navier-Stokes equations [23]. The model equations that represent the movement of drug particles flowing through the nasal cavity or lungs include a system of equations for a continuous medium, and an equation for a dispersed phase, which describes the movement of particles and includes the sum of all forces that affect their movement. The system of equations for a continuous medium (1) includes the continuity equation, the momentum conservation equation, and the energy conservation equation [24]: where α and α p are the proportion of the continuous medium and the dispersed phase; ρ is the density of the continuous medium, kg/m 3 ; v → is the velocity vector of the continuous medium, m/s; v → p is the particle velocity vector, m/s; p is the static pressure, Pa; τ is the viscous stress tensor; ρg → is the gravity force, kg/m 2 ·s 2 ; R → sl is the force that characterizes the effect of the dispersed phase on the continuous medium, kg/m 2 ·s 2 ; f → p is the force characterizing the resistance of the particle, N; V p is the particle volume, m 3 ; K sl is the momentum exchange between phases, kg/(m 3 ·s); µ is the dynamic viscosity, Pa·s; C P is the specific heat capacity of the mixture, J/kg·K; T is the temperature of the mixture, K; λ is the thermal conductivity of the mixture, W/m·K; δ k I is the unit tensor; and i is the particle index.
To describe the motion of the dispersed phase, Equation (2) [25] is used: In this case, the interaction between a particle and a continuous medium is described by the following equations: where m p is the particle mass, kg; g → is the free fall acceleration, m/s 2 ; f → p, c is the force of interaction between a continuous medium and a particle, N; f → p, p is the force of interaction between particles, N; f → p, w is the force of interaction between the particle and the wall, N; f → ∇p is the pressure gradient force acting on a particle, N; ∇p is the pressure gradient, Pa; f → ∇τ is the viscous force, N; Num is the number of particles, i; and j is the particle indices.
In addition to the equations described above, in a number of works [24,25], interactions between colliding particles are taken into account; for this, the stiffness coefficient is calculated, on the basis of which the forces of interaction between particles are calculated.
Depending on the region of the respiratory system under consideration, a change in the nature of the movement of the air flow is observed. In [23], the authors demonstrated the presence of turbulent and transient flows in the extrathoracic region. The presence of turbulent air flow in the nasal and oral cavities can have a significant impact on the nature of the movement and the deposition zone of drug particles; however, a number of works [26][27][28] have shown that there is no turbulent flow in the tracheobronchial and alveolar regions.
The described features of the air flow in various areas of the respiratory system necessitate the correct choice of the method for calculating turbulent flows in the study of motion and drug deposition zones in silico using computational fluid dynamics.
At the moment, there are a large number of methods for calculating turbulent flows based on the representation of a turbulent flow as a spectrum of eddies of various scales: the method of direct numerical simulation (DNS), the method of large eddies (LES) and the method of averaging the Navier-Stokes equations according to Reynolds (RANS). One of the priority areas for calculating turbulent flows is numerical simulation based on the construction of difference calculation methods.
A significant place in modern research is occupied by the DNS calculation method, which allows for solving the complete non-stationary, three-dimensional Navier-Stokes equations without using special turbulence models, i.e., models k-ε, and k-ω. The complexity of the DNS method is primarily due to the fact that unsteady turbulent flows are characterized by a wide range of spatial and temporal scales. Turbulence scales are distributed over a range of linear dimensions-from the largest eddies interacting with an average flow, the time scale of which is comparable to the scale of large eddies, to the smallest eddies, called the Kolmogorov scale, dissipating under the action of viscous forces. To perform calculations using the DNS method, the computational grid must be fine enough to resolve the smallest eddies, the dimensions of which are of the order of the Kolmogorov length scale: where η k is the Kolmogorov scale, m/s; ν is the kinematic viscosity, m 2 /s; and ε is the energy dissipation, m 2 /s 3 . The DNS method is applicable only for limited Reynolds numbers (Re), because the number of N DNS grid nodes required for calculations is proportional to the following ratio: Using relation (8), one can estimate the number of N DNS grid nodes. At Reynolds numbers of the order of 10 4 , the number of nodes of the computational grid will be 1.4 × 10 7 , and with the Reynolds number equal to 10 5 , the number of nodes of the computational grid will be 2. 49 10 9 . Therefore, calculations by the DNS method require high-performance computing equipment and an efficient numerical method that makes it possible to obtain reliable numerical results. Such requirements limit the widespread use of the DNS method in describing the movement of drug particles in the respiratory system. In [29], the emergence of a turbulent vortex in the oral cavity was studied; the DNS method was used to calculate turbulent flows. The data obtained indicate the presence of high turbulence in the oral cavity and larynx.
The presence of turbulent vortices in the oral cavity and larynx was proved in [30,31] using the RANS method based on the replacement of changing flow characteristics (i.e., velocity, pressure, and density) by sums of averaged and pulsation components, using semi-empirical turbulence models. The approaches proposed in the papers have demonstrated a significant reduction in computing power by reducing the requirements for the discretization of the computational domain and the time step. In addition, the LES method [32,33] can also be used to calculate turbulent eddies.
The RANS (Reynolds-Averaged Navier-Stokes) calculation method is commonly used to study the airflow dynamics in the respiratory system due to its adequate accuracy and low computational requirements; however, in order to solve the system of Reynoldsaveraged Navier-Stokes equations, additional empirical relationships are needed, specifically, rheological models of states. These relationships are based on hypotheses related to the turbulent viscosity, mixing paths, and turbulent flows.
It is worth noting that when modeling the movement and deposition of particles or droplets of drugs, it is essential to consider not only the methods for calculating turbulent flows but also models that describe the motion and interaction of particles with each other, the continuous phase, and the walls of the respiratory system or nasal cavity.

Methods for Describing the Movement of Particles or Droplets of Drugs in the Respiratory System
The Euler and Lagrange models can be used to describe the movement of particles or droplets of drugs in the inhalation and intranasal routes of administration; however, when implementing the Euler model, drug particles or drops are considered as one continuous phase. In turn, the behavior of the flow is predicted based on the resolution of the laws of conservation of mass and momentum.
The Lagrange model allows us to consider drug particles or drops as a separate dispersed phase in a continuous gaseous medium. Particle motion is predicted by integrating Newton's second law with respect to time. The Lagrange model makes it possible to study the dynamics of the movement of particles or drops and the zone of their deposition.
Depending on the consideration of the interaction of particles of the dispersed phase with each other, this approach can be divided into the discrete particle method (DPM), which does not consider the particle-particle, or the particle-wall interaction, and the discrete element method (DEM), which takes into account these interactions. When using the discrete particle method (DPM), it is assumed that the volume fraction of the dispersed phase in a continuous medium is insignificant; therefore, it does not affect the nature of the air flow, and particle-particle interactions are not taken into account [34]. In addition, it is assumed that the precipitation of particles or droplets of a drug occurs upon contact with the walls of the organs of the respiratory system. In most studies, in order to carry out mathematical modeling and achieve the required accuracy of calculations, the method of discrete particles has been used when considering the forces of resistance and gravity [35][36][37]. The authors of [38] evaluated the effect of gravitational forces on the zones of deposition of drug particles. It was established that in the absence of gravity, the total deposition of particles in the tracheobronchial region of the respiratory system was reduced by 10%. This was due to an increase in the proportion of precipitated drug particles in the area of administration (i.e., the oral cavity).
In turn, at a high concentration of particles in the dispersed phase, about 10 4 or 5·10 3 particles/mL of gas for aerosols and 1.2·10 9 particles per dose for inhalation powders [39,40], it is necessary to take into account the interaction between particles, between particles and the flow continuous medium, and between particles and the walls of the respiratory system organs. The presented interactions can be taken into account when using computational fluid dynamics in conjunction with the discrete element method [41]. In this case, particle-particle, particle-wall, and particle-continuous medium interactions are taken into account. This approach is widely used in the development and evaluation of the effectiveness of new inhalation and intranasal drug delivery systems; however, there are limitations to using this approach for the mathematical modeling of particle movement and deposition zones.
The authors of [42] carried out mathematical modeling of particle motion using the method of computational fluid dynamics in combination with the method of discrete elements in the Weibel light model. The data obtained indicated that the initial parameters of the introduction of particles have a significant impact on their movement in the respiratory system. In [43], the authors compared two approaches to modeling particle motion (i.e., the discrete particle method and the discrete element method) in the extrathoracic region of the respiratory system. When using the discrete element method, a complex trajectory of particles in the upper part of the trachea was revealed, due to the strong interaction between the particles. Therefore, when employing the discrete element method to simulate particle movement, a notable deposition of particles in the extrathoracic region is observed. These findings are consistent with the studies presented in [44].

Influence of Flow Characteristics on the Deposition Zones and Movement of Medicinal Substances
Airflow characteristics play a crucial role in determining flow patterns and drug deposition zones. Currently, substantial advancements have been made in the research focused on investigating the impact of airflow parameters on particle deposition zones [45][46][47].

1.
Effect of flow rate.
Increasing the velocity results in turbulence and the inertial settling of drug particles, which causes deposition in the extrathoracic region of the respiratory system. In [48], the authors studied the transition from a laminar to turbulent flow in the respiratory system. It was found that an increase in the flow velocity leads to the formation of turbulent vortices in the oral cavity. In turn, fluctuations in the flow velocity have a significant effect on the movement of particles and their deposition. The authors of [49] analyzed the effect of the inhalation velocity of drug particles, equal to 30 and 60 L/min (3.6 and 7.2 m/s, respectively, with the cross-sectional area of the inlet opening equal to 0.00014 m 2 ), corresponding to pharmacopeia data, on their deposition zones using computational fluid dynamics and a lung model built on the basis of computed tomography. It was found in the work that with an increase in the flow rate, a significant settling of particles in the oral cavity was observed. However, the authors of [50] found that the flow velocity does not affect the deposition zones below the extrathoracic region.
In turn, the authors of [51] carried out CFD modeling of the flow with an inlet velocity of 28.3 and 60 L/min in a cascade Andersen impactor divided into 7 stages, each of which had a characteristic hole diameter used to study the aerodynamic distribution of fine particles. Based on the simulation results, it was found that at the impactor stages 0 and 1 (i.e., a hole diameter of 2.55-1.89 mm), the flow was divided in directions due to the central hole, the formation of small eddy currents around the impactor outlet holes was observed, which were most noticeable at stages 2 and 3 (i.e., a hole diameter of 0.91-0.71 mm) and they become less pronounced in steps 4 to 7 (i.e., a hole diameter of 0.53-0.25 mm). Table 1 shows the characteristic particle diameters deposited on the respective impactor stages at different flow rates. The authors of [51] showed that at a higher inlet flow rate of 60 L/min, the characteristic particle size deposited in the first three stages of the impactor is less than at a rate of 28.3 L/min.
To study the effect of flow velocity on particle deposition zones, the authors of [52] used casts of the human upper respiratory tract with a subsequent digitization of the periphery of the airway contour (Figure 4a). The authors of [52] measured the crosssectional area and perimeter of the oral cavity, pharynx, and larynx by cutting the cast into 0.3 cm sections (Figure 4b). The authors of [51] showed that at a higher inlet flow rate of 60 L/min, the characteristic particle size deposited in the first three stages of the impactor is less than at a rate of 28.3 L/min.
To study the effect of flow velocity on particle deposition zones, the authors of [52] used casts of the human upper respiratory tract with a subsequent digitization of the periphery of the airway contour (Figure 4a). The authors of [52] measured the cross-sectional area and perimeter of the oral cavity, pharynx, and larynx by cutting the cast into 0.3 cm sections (Figure 4b). Using an airway model, the authors of [53] studied the effect of various flow rates, namely, 15, 30 and 60 L/min, on the deposition efficiency of particles of various sizes in the respiratory system. The results of the studies of the authors of [53] are presented in Figure 5a,b. Using an airway model, the authors of [53] studied the effect of various flow rates, namely, 15, 30 and 60 L/min, on the deposition efficiency of particles of various sizes in the respiratory system. The results of the studies of the authors of [53] are presented in Figure 5a,b. The results showed that the particle settling increased with an increasing flow rate. At the highest flow rate (60 L/min), more than 90% of the particles, with a characteristic size >20 µm, were deposited in the oral cavity (Figure 5a). At a flow rate of 30 L/min for particles <10 µm, a more uniform settling pattern was observed in three regions of the respiratory system: the oral cavity, pharyngeal and laryngeal regions, and tracheobronchial (TB), while the settling efficiency did not exceed 10%; however, for particles >10 µm, the deposition efficiency in the oral cavity was higher (Figure 5b).
At the same time, the authors of work [54] note that in order to perform correct calculations, the number of cells of the computational grid should be at least 480,000 for flow rates of 15 and 30 L/min, and 640,000 for a flow rate of 60 L/min.
To study the particle settling zone, different respiratory modes are compared, which affect the flow rate. The authors of works [55][56][57][58] carried out mathematical modeling us- The results showed that the particle settling increased with an increasing flow rate. At the highest flow rate (60 L/min), more than 90% of the particles, with a characteristic size >20 µm, were deposited in the oral cavity (Figure 5a). At a flow rate of 30 L/min for particles <10 µm, a more uniform settling pattern was observed in three regions of the respiratory system: the oral cavity, pharyngeal and laryngeal regions, and tracheobronchial (TB), while the settling efficiency did not exceed 10%; however, for particles >10 µm, the deposition efficiency in the oral cavity was higher (Figure 5b). At the same time, the authors of work [54] note that in order to perform correct calculations, the number of cells of the computational grid should be at least 480,000 for flow rates of 15 and 30 L/min, and 640,000 for a flow rate of 60 L/min.
To study the particle settling zone, different respiratory modes are compared, which affect the flow rate. The authors of works [55][56][57][58] carried out mathematical modeling using computational fluid dynamics to study the influence of a non-stationary and stationary inspiratory velocity flow for a simplified model of the oral cavity and larynx ( Figure 6).
The results showed that the particle settling increased with an increasing flow rate. At the highest flow rate (60 L/min), more than 90% of the particles, with a characteristic size >20 µm, were deposited in the oral cavity (Figure 5a). At a flow rate of 30 L/min for particles <10 µm, a more uniform settling pattern was observed in three regions of the respiratory system: the oral cavity, pharyngeal and laryngeal regions, and tracheobronchial (TB), while the settling efficiency did not exceed 10%; however, for particles >10 µm, the deposition efficiency in the oral cavity was higher (Figure 5b).
At the same time, the authors of work [54] note that in order to perform correct calculations, the number of cells of the computational grid should be at least 480,000 for flow rates of 15 and 30 L/min, and 640,000 for a flow rate of 60 L/min.

Influence of unsteady flow.
To study the particle settling zone, different respiratory modes are compared, which affect the flow rate. The authors of works [55][56][57][58] carried out mathematical modeling using computational fluid dynamics to study the influence of a non-stationary and stationary inspiratory velocity flow for a simplified model of the oral cavity and larynx ( Figure  6). The results presented in the work indicate the presence of vortex flows in the oral cavity, larynx and trachea region both in a stationary and non-stationary air flow; however, with an unsteady flow, the vortex motion takes on a more complex character. The results presented in the work indicate the presence of vortex flows in the oral cavity, larynx and trachea region both in a stationary and non-stationary air flow; however, with an unsteady flow, the vortex motion takes on a more complex character.
In addition, a study conducted by researchers in [57] demonstrated that the presence of an unsteady flow, characterized by time-dependent velocity and pressure fields, leads to the increased deposition of medicinal substances. Moreover, when the airflow exhibits a cyclic nature corresponding to the inhalation and exhalation cycles, it was observed that turbulent flow develops during the acceleration phase of inspiration, followed by a transition to a laminar flow as the flow decelerates [58].

3.
Influence of the physiological model of breathing.
In [59] it was demonstrated that the nature of breathing, especially during inhalation, has a significant effect on the movement and deposition of medicinal substances. The effect of the respiration cycle on the deposition and movement of medicinal substances was also studied by the authors of [60]. It was concluded that with each breathing cycle, the concentration of drugs deposited in the organs of the respiratory system increases, and subsequent breathing cycles lead to the further movement of particles to the lower parts of the respiratory system. When considering the delivery of drugs based on nanoparticles, it was found that, towards the end of inhalation, the efficiency of particle deposition increases with a decreasing flow rate due to an increase in the residence time of the particles [61].
For some dosing devices, it is necessary to hold the breath in order to precipitate suspended particles of the medicinal substance in the organs of the respiratory system. The authors of [62], using computational fluid dynamics, studied the movement and deposition of drug particles in the respiratory tract from the oral cavity to the bronchi of the 7th order.
The data obtained indicate the presence of turbulence during breath holding. It was found in the work that for particles with an average size of 5 µm, holding the breath makes it possible to reduce the probability of deposition of particles in the oral cavity with an increase in deposition in the bronchi. Thus, for a number of drug delivery systems, the deposition of which is necessary in the lower respiratory tract (Figure 2), it is necessary to hold the breath.

Computational Fluid Dynamics for Studying the Motion and Deposition of Medicinal Substances
Studying the movement and deposition of drug particles using computational fluid dynamics and in vivo experimental studies is an important tool in the development of new innovative inhalation and intranasal drug delivery systems. For the first time, the use of computational fluid dynamics for the development of an inhalation delivery system was presented by the authors of [63]. Various dosing device designs and continuum compositions were investigated to achieve the required area of drug particle deposition in the respiratory system.
The data obtained indicate that the reduction in the metering device valve outlets and the presence of a spacer allow the delivery of drugs to the bronchi without deposition in the oral cavity.
In [64], the possibility of a targeted delivery of inhaled drugs using magnetic targeting was studied. The authors conducted a computational experiment using computational fluid dynamics with the method of calculating turbulent flows with averaging the Navier-Stokes equations according to Reynolds (Figure 7). The authors of [64] showed that the local and overall efficiency of particle deposition can be significantly increased by activating the strength of the magnetic field. The data obtained indicate the possible further development of the targeted delivery of drugs through magnetic targeting for the treatment of diseases of the respiratory system.
Based on the considered approaches to modeling the movement and deposition of drug particles in various lung segments, the main conclusions can be drawn. First, depending on the region of the respiratory system under consideration, a change in the nature of the movement of the air flow is observed; therefore, to calculate the movement and deposition of drug particles, it is necessary to take into account the occurrence of turbulent flow-the most appropriate and widely used calculation method, in the opinion of the authors, is RANS, since the hypotheses of turbulent viscosity, mixing paths and turbulent flows, and the computing power requirements are lower than for the DNS method. Secondly, it was shown that at a high concentration of particles or drops in the dispersed phase, it is necessary to consider not only methods for calculating the turbulent flows, but also models of the movement and interaction of particles or drops with each other, the continuous phase and the walls of the respiratory system. Thirdly, the considered influence of the air flow parameters on the particle deposition zones made it possible to determine that an increase in the inlet flow rate contributes to the deposition of a significant The authors of [64] showed that the local and overall efficiency of particle deposition can be significantly increased by activating the strength of the magnetic field. The data obtained indicate the possible further development of the targeted delivery of drugs through magnetic targeting for the treatment of diseases of the respiratory system.
Based on the considered approaches to modeling the movement and deposition of drug particles in various lung segments, the main conclusions can be drawn. First, depending on the region of the respiratory system under consideration, a change in the nature of the movement of the air flow is observed; therefore, to calculate the movement and deposition of drug particles, it is necessary to take into account the occurrence of turbulent flow-the most appropriate and widely used calculation method, in the opinion of the authors, is RANS, since the hypotheses of turbulent viscosity, mixing paths and turbulent flows, and the computing power requirements are lower than for the DNS method. Secondly, it was shown that at a high concentration of particles or drops in the dispersed phase, it is necessary to consider not only methods for calculating the turbulent flows, but also models of the movement and interaction of particles or drops with each other, the continuous phase and the walls of the respiratory system. Thirdly, the considered influence of the air flow parameters on the particle deposition zones made it possible to determine that an increase in the inlet flow rate contributes to the deposition of a significant amount of particles in the oral cavity; however, by choosing the right dosing device it is possible to achieve the deposition of drug particles directly in the bronchi without deposition in the oral cavity. In addition, the presence of breathing cycles (i.e., inhalation and exhalation) also affects the efficiency of the delivery of drug particles to the lungs. A decrease in the flow rate is observed towards the end of inspiration, which leads to an increase in the residence time of the particles; therefore, the efficiency of their deposition increases. The authors of this paper believe that the considered approaches to mathematical modeling of the trajectory of movement and deposition of particles or drops of a drug will contribute to the accelerated development of new inhaled drugs needed to treat common lung diseases.

Computational Fluid Dynamics for Studying the Movement and Deposition of Medicinal Substances in the Nasal Cavity
Currently, the development of intranasal drug delivery systems for the treatment of various neurological disorders is relevant and promising [65]. The onset of the desired therapeutic effect is possible when it is reached (Figure 8). Existing approaches, however, do not allow for achieving a high concentration of particles in the required area. This is due to both the complex structure of the nasal cavity, in which inhaled particles are filtered long before reaching the olfactory region, and the complexity of controlling the trajectory of drug particles. Various approaches to intranasal drug administration have been studied to improve delivery efficiency to the olfactory region.
In [66], using computational fluid dynamics, a numerical simulation of the flow trajectory was carried out when a drug was injected into a certain area of the nasal cavity. The obtained results indicate that the introduction of drugs into the region of the anterior wall of the nasal cavity increases the deposition coefficient of particles in the olfactory region compared to the traditional method of administration.
The authors of [67], using the method of computational fluid dynamics in an anatomically-accurate model of the nose, studied the efficiency of delivering a group of charged aerosols to the osteomeatal complex (OMC) when exposed to an electric field. The work considered both monodisperse and polydisperse aerosol profiles. Figure 9 schematically shows the effect of the electric field on the deposition zones of drug particles during the intranasal route of administration [67]. Existing approaches, however, do not allow for achieving a high concentration of particles in the required area. This is due to both the complex structure of the nasal cavity, in which inhaled particles are filtered long before reaching the olfactory region, and the complexity of controlling the trajectory of drug particles. Various approaches to intranasal drug administration have been studied to improve delivery efficiency to the olfactory region.
In [66], using computational fluid dynamics, a numerical simulation of the flow trajectory was carried out when a drug was injected into a certain area of the nasal cavity. The obtained results indicate that the introduction of drugs into the region of the anterior wall of the nasal cavity increases the deposition coefficient of particles in the olfactory region compared to the traditional method of administration.
The authors of [67], using the method of computational fluid dynamics in an anatomicallyaccurate model of the nose, studied the efficiency of delivering a group of charged aerosols to the osteomeatal complex (OMC) when exposed to an electric field. The work considered both monodisperse and polydisperse aerosol profiles. Figure 9 schematically shows the effect of the electric field on the deposition zones of drug particles during the intranasal route of administration [67]. jectory was carried out when a drug was injected into a certain area of the nasal cavity. The obtained results indicate that the introduction of drugs into the region of the anterior wall of the nasal cavity increases the deposition coefficient of particles in the olfactory region compared to the traditional method of administration.
The authors of [67], using the method of computational fluid dynamics in an anatomically-accurate model of the nose, studied the efficiency of delivering a group of charged aerosols to the osteomeatal complex (OMC) when exposed to an electric field. The work considered both monodisperse and polydisperse aerosol profiles. Figure 9 schematically shows the effect of the electric field on the deposition zones of drug particles during the intranasal route of administration [67]. Mathematical modeling made it possible to determine that the synthesis of electrical guidance and the point release of the drug allows for achieving a focused deposition of aerosol droplets with a significantly increased dosage in the OMC. Table 2 shows the mass fractions of deposited monodisperse particles in the OMC region. Mathematical modeling made it possible to determine that the synthesis of electrical guidance and the point release of the drug allows for achieving a focused deposition of aerosol droplets with a significantly increased dosage in the OMC. Table 2 shows the mass fractions of deposited monodisperse particles in the OMC region. Meanwhile, Table 3 shows the mass fractions of deposited polydisperse particles in the OMC region. The results of the performed mathematical modeling (Tables 2 and 3) showed the delivery efficiency in the OMC to be 51.6% for monodisperse aerosols and 34.4% for polydisperse aerosols. These differences indicate that the aerosol profile has a significant impact on intranasal drug administration.
The considered approaches to mathematical modeling using computational fluid dynamics can be useful in the development of new intranasal drug delivery systems needed to treat a wide range of diseases from rhinosinusitis to neurological and mental disorders. In the authors' opinions, however, mathematical modeling of intranasal drug delivery using new approaches (e.g., exposure to electric and electromagnetic fields), which increase the delivery efficiency and drug concentration in the target area of the nasal cavity, should be carefully tested by conducting a series of in vitro experiments.
In the development of innovative systems for inhalation and intranasal drug delivery, studies aimed at multiscale modeling are becoming increasingly important. In order to evaluate the effectiveness of the developed delivery systems, computational fluid dynamics is used to study the areas of deposition of particles or drops of drugs in the respiratory system, using various dissolution models, namely, the kinetics of dissolution, and the release of drugs in a specific area of the respiratory system are predicted, while the results obtained are used to build pharmacokinetic models to describe the distribution of an active pharmaceutical ingredient throughout the body depending on the area of deposition [68].

Description of the Kinetics of the Release of the Active Substance from the Dosage Form
After the precipitation of particles or droplets of drugs in the respiratory system, the active substance is released from the dosage form. The release is a complex process that sequentially includes the following stages: the interaction of the active substance with the environment in which the precipitation of the drug particles occurred; dissolution of the active substance; and the movement of the dissolved active substance within the dosage form and, then, in the human body.
Thus, the task of modeling the drug release process includes the following steps: 1.
Prediction of drug particle deposition zones in the respiratory system.

2.
Modeling the internal heterogeneous structure of a dosage form with an embedded active substance.

3.
Modeling the process of dissolution of the active substance.

4.
Simulation of the active substance's transportation within the dosage form, considering the medium's heterogeneity, and subsequent transportation outside the dosage form.
Mathematical models for the dissolution process are based on the solution of differential equations. Their fundamental principle is to examine the entire system and track the changes in the concentration of the dissolved active substance over time.
However, a significant drawback of these models is their inaccurate description of the dissolution process for multicomponent solids and bodies with complex internal structures. This limitation stems from various factors and the stochastic nature of the multicomponent dosage form's structure. Therefore, the development of models that consider the complex and heterogeneous internal structure of the dosage form remains an urgent task in the pursuit of a comprehensive multiscale model.

Modeling the Internal Structure of a Dosage form with an Embedded Active Substance and Release from It
Modeling the structure of a dosage form can be considered at several hierarchical levels: the molecular level, micro level and macro level, each of which has its own methods and approaches, as well as limitations.
As part of the development of a single multiscale model, the most promising are approaches that consider the system at the micro level, since they are still sufficiently detailed to take into account the geometry of the structure of the dosage form, in contrast to the macro level, but they do not require such large computing power as approaches related to the molecular level.
At the micro level, the main modeling methods are the method of dissipative particle dynamics, the method of lattice Boltzmann equations, the method of Brownian dynamics, and the method of a self-consistent field. The system in this case is represented by describing the field or conditional microscopic particles.
In the case of heterogeneous nanoporous structures, such as aerogel dosage forms, it is often necessary to apply special modeling approaches. This is due to the fact that the shape of such structures is stochastic in nature and cannot be described by geometric objects [76].
One of the promising approaches for describing nanoporous systems at the microlevel is the cellular automaton (CA) approach.
The essence of the CA approach is to divide the system into identical local areas (cells), while each cell at any given time can only be in one of a finite, predetermined set of states. In this case, time is also divided into equal discrete steps, and the change in the state of the cell at each time step is determined in accordance with the specified transition rules. This approach allows us to move away from a complex system of equations, in which many factors must be taken into account, to local rules of cell behavior that do not depend on each other. This makes it possible to significantly simplify the mathematical apparatus in the calculation of complex systems while maintaining the required accuracy. This makes the CA approach a suitable tool for describing complex heterogeneous systems such as nanoporous dosage forms [77], various kinds of coated or uncoated polymer particles, and polymer particles embedded in an active pharmaceutical ingredient (API). Figure 10 shows examples of modeling various structures of drug particles based on the cellular automaton approach. It should be noted that the CA approach can also be used in modeling the dissolution of an active substance and is a promising alternative to the traditional approach, which consists in developing and solving differential equations. In particular, the CA approach has been successfully used in modeling the dissolution of the active substance from solid dosage forms. In [77], a model for the dissolution of three-dimensional solids was proposed based on the cellular automaton approach.
In addition, one of the advantages of the CA approach is the ability to combine different CA models. In [7], a CA model was implemented that described the release of melatonin in a model nasal cavity during the deposition of aerogel particles based on chitosan. Thus, it was possible to develop a unified CA model of the structure of the dosage form and the process of dissolution and transportation of the active substance in it.

Existing Approaches for Modeling the Movement of the Active Substance after Its Release
Simulation of the transport of an active substance in solution is the problem of modeling the hydrodynamics of a multiphase multicomponent system.
Currently, the most common methods for modeling hydrodynamics rely on solving the system of Navier-Stokes differential equations; however, these methods have limited applicability when dealing with multi-phase and multi-component systems, and they are also computationally demanding.
The most suitable method for developing a multiscale model, which allows for simplified mathematical modeling and reduced computational costs, is the lattice Boltzmann method based on lattice Boltzmann equations [78].
The main idea behind the lattice Boltzmann method is to divide the system into identical cells. Each cell represents a volume of the simulated medium containing fictitious fluid particles. Fluid particles can only move between adjacent cells in each time step. The It should be noted that the CA approach can also be used in modeling the dissolution of an active substance and is a promising alternative to the traditional approach, which consists in developing and solving differential equations. In particular, the CA approach has been successfully used in modeling the dissolution of the active substance from solid dosage forms. In [77], a model for the dissolution of three-dimensional solids was proposed based on the cellular automaton approach.
In addition, one of the advantages of the CA approach is the ability to combine different CA models. In [7], a CA model was implemented that described the release of melatonin in a model nasal cavity during the deposition of aerogel particles based on chitosan. Thus, it was possible to develop a unified CA model of the structure of the dosage form and the process of dissolution and transportation of the active substance in it.

Existing Approaches for Modeling the Movement of the Active Substance after Its Release
Simulation of the transport of an active substance in solution is the problem of modeling the hydrodynamics of a multiphase multicomponent system. Currently, the most common methods for modeling hydrodynamics rely on solving the system of Navier-Stokes differential equations; however, these methods have limited applicability when dealing with multi-phase and multi-component systems, and they are also computationally demanding.
The most suitable method for developing a multiscale model, which allows for simplified mathematical modeling and reduced computational costs, is the lattice Boltzmann method based on lattice Boltzmann equations [78].
The main idea behind the lattice Boltzmann method is to divide the system into identical cells. Each cell represents a volume of the simulated medium containing fictitious fluid particles. Fluid particles can only move between adjacent cells in each time step. The direction of particle movement in the next time step is strictly defined. The calculation of each time step occurs in two phases: the flow propagation phase and the collision phase.
In [7], the movement of fluid particles to neighboring cells was considered in detail in accordance with their direction (Figure 11). During the first phase, particles move from one cell to another in accordance with their directions. The number of possible directions of particle motion depends on the number of dimensions of the system (i.e., a two-dimensional or three-dimensional lattice), the required accuracy, and the computational power [79].
The flow propagation phase is followed by the collision phase. At this step, the collision of particles with each other inside the cell is calculated, thereby changing the distribution of particles along the directions of movement. How a collision is handled is determined by the collision operator; therefore, the distribution after collisions will be different.
The method of lattice Boltzmann equations has shown itself well in modeling multicomponent systems with complex geometry [80]. In addition, it is also used in medical research, in particular in [81], where a variant of the model is proposed, which is suitable for modeling the blood flow in vessels.
This conducted literature review has enabled the identification of the most pressing tasks and directions in the field of mathematical modeling of the respiratory system. The combination of computational fluid dynamics approaches to calculate the trajectory of particle movement and deposition in the respiratory system, along with the cellular automaton (CA) approach to model the kinetics of drug dissolution and release from the dosage form, followed by the calculation of drug redistribution throughout the body via blood flow based on the deposition area [82][83][84][85][86][87][88][89], will facilitate the development of a unified multiscale model.
Based on the collected scientific and technical information presented in this review, at the moment it is possible to propose the following type of a single multi-scale model, as a new approach for the accelerated development of inhaled and intranasal drugs necessary for the treatment of socially significant diseases- Figure 12. During the first phase, particles move from one cell to another in accordance with their directions. The number of possible directions of particle motion depends on the number of dimensions of the system (i.e., a two-dimensional or three-dimensional lattice), the required accuracy, and the computational power [79].
The flow propagation phase is followed by the collision phase. At this step, the collision of particles with each other inside the cell is calculated, thereby changing the distribution of particles along the directions of movement. How a collision is handled is determined by the collision operator; therefore, the distribution after collisions will be different.
The method of lattice Boltzmann equations has shown itself well in modeling multicomponent systems with complex geometry [80]. In addition, it is also used in medical research, in particular in [81], where a variant of the model is proposed, which is suitable for modeling the blood flow in vessels.
This conducted literature review has enabled the identification of the most pressing tasks and directions in the field of mathematical modeling of the respiratory system. The combination of computational fluid dynamics approaches to calculate the trajectory of particle movement and deposition in the respiratory system, along with the cellular automaton (CA) approach to model the kinetics of drug dissolution and release from the dosage form, followed by the calculation of drug redistribution throughout the body via blood flow based on the deposition area [82][83][84][85][86][87][88][89], will facilitate the development of a unified multiscale model.
Based on the collected scientific and technical information presented in this review, at the moment it is possible to propose the following type of a single multi-scale model, as a new approach for the accelerated development of inhaled and intranasal drugs necessary for the treatment of socially significant diseases- Figure 12. The data collected from the literature analysis will be utilized by the authors of this study in the development of a multiscale model. The objective of this model is to expedite the development and introduction of new inhaled and intranasal drugs required for the treatment of socially significant diseases into the market. The data collected from the literature analysis will be utilized by the authors of this study in the development of a multiscale model. The objective of this model is to expedite the development and introduction of new inhaled and intranasal drugs required for the treatment of socially significant diseases into the market.

Conclusions
A review of the scientific and technical literature in the field of mathematical modeling showed that, at the moment, there is considerable experience in using models of different levels to describe certain classes of forecasting problems: flow hydrodynamics during spraying and the further movement of particles or drops, as well as the area of their deposition; the kinetics of dissolution and release of the drug substance; and the adsorption, redistribution and excretion of a medicinal substance in the organs and tissues of the body. The research results presented in this article give a broad idea of the problems that arise in the mathematical modeling of anatomically-similar organs of the respiratory system, and they draw attention to various factors that affect the relationship between the quality of the computational grid and the computation time. The presented studies in the field of CFD modeling allow for choosing the most optimal calculation methods suitable for the organs of the respiratory system. In addition, the results of the studies draw attention to promising methods for the targeted delivery of drugs to the organs of the human respiratory system. With the help of mathematical modeling tools, it is shown how a particular method affects the efficiency of drug delivery in specific zones/segments of the respiratory system. The presented results in the field of cellular automaton models show that this mathematical modeling tool is promising in the field of drug development, due to simplified mathematical modeling and a reduced computational time. Combining models of different levels also gives researchers a comprehensive tool for preliminary reverse engineering in development. Nevertheless, despite significant advances in the field of mathematical modeling, the parametrization of models and verification of their adequacy remain a current problem, since complex models require a significant amount of experimental data at the stages of parameterization and verification. Modern medicine today does not have a non-invasive, commercially available and accurate tool for obtaining complete information about the transformation of the delivery system and drug substances inside the body. The developers of mathematical models have to work with disparate blocks of data coming from different sources, both at the stage of preclinical trials and from the stages of clinical trials. A significant part of the parameters of complex models is currently determined by optimization methods, and they do not guarantee a single correct solution, and not from experimental procedures. The problem of the lack of complete information hinders the widespread use of mathematical models as one of the main reverse engineering tools in the pharmaceutical development process. The collected literature review should be the first step towards the development of a single multiscale model aimed at the development of new drugs needed for the treatment of socially significant diseases.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article.