All issues
- 2025 Vol. 17
- 2024 Vol. 16
- 2023 Vol. 15
- 2022 Vol. 14
- 2021 Vol. 13
- 2020 Vol. 12
- 2019 Vol. 11
- 2018 Vol. 10
- 2017 Vol. 9
- 2016 Vol. 8
- 2015 Vol. 7
- 2014 Vol. 6
- 2013 Vol. 5
- 2012 Vol. 4
- 2011 Vol. 3
- 2010 Vol. 2
- 2009 Vol. 1
-
Optimization of proton therapy with radiosensitizing nanoparticles and antiangiogenic therapy via mathematical modeling
Computer Research and Modeling, 2025, v. 17, no. 4, pp. 697-715Optimization of antitumor radiotherapy represents an urgent issue, as approximately half of the patients diagnosed with cancer undergo radiotherapy during their treatment. Proton therapy is potentially more efficient than traditional X-ray radiotherapy due to fundamental differences in physics of dose deposition, leading to better targeting of tumors and less collateral damage to healthy tissue. There is increasing interest in the use of non-radioactive radiosensitizing tumor-specific nanoparticles the use of which can boost the performance of proton therapy. Such nanoparticles are small volumes of a sensitizer, such as boron-10 or various metal oxides, enclosed in a polymer layer containing tumor-specific antibodies, which allows for their targeted delivery to malignant cells. Furthermore, a combination of proton therapy with antiangiogenic therapy that normalizes tumor-associated microvasculature may yield further synergistic increase in overall treatment efficacy.
We have developed a spatially distributed mathematical model simulating the growth of a non-invasive tumor undergoing treatment by fractionated proton therapy with nanosensitizers and antiangiogenic therapy. The modeling results suggest that the most effective way to combine these treatment modalities should strongly depend on the tumor cells’ proliferation rate and their intrinsic radiosensitivity. Namely, a combination of antiangiogenic therapy with proton therapy, regardless of whether radiosensitizing nanoparticles are used, benefits treatment efficacy of rapidly growing tumors as well as radioresistant tumors with moderate growth rate. In these cases, administration of proton therapy simultaneously with antiangiogenic drugs after the initial single injection of nanosensitizers is the most effective option among those analyzed. Conversely, for slowly growing tumors, maximization of the number of nanosensitizer injections without antiangiogenic therapy proves to be a more efficient option, with enhancement in treatment efficacy growing with the increase of tumor radiosensitivity. However, the results also show that the overall efficacy of proton therapy is likely to increase only modestly with the addition of nanosensitizers and antiangiogenic drugs.
Keywords: mathematical oncology, numerical optimization. -
A minimal model of density-dependent population dynamics incorporating sex structure: simulation and application
Computer Research and Modeling, 2025, v. 17, no. 5, pp. 941-961This study proposes and analyzes a discrete-time mathematical model of population dynamics with seasonal reproduction, taking into account the density-dependent regulation and sex structure. In the model, population birth rate depends on the number of females, while density is regulated through juvenile survival, which decreases exponentially with increasing total population size. Analytical and numerical investigations of the model demonstrate that when more than half of both females and males survive, the population exhibits stable dynamics even at relatively high birth rates. Oscillations arise when the limitation of female survival exceeds that of male survival. Increasing the intensity of male survival limitation can stabilize population dynamics, an effect particularly evident when the proportion of female offspring is low. Depending on parameter values, the model exhibits stable, periodic, or irregular dynamics, including multistability, where changes in current population size driven by external factors can shift the system between coexisting dynamic modes. To apply the model to real populations, we propose an approach for estimating demographic parameters based on total abundance data. The key idea is to reduce the two-component discrete model with sex structure to a delay equation dependent only on total population size. In this formulation, the initial sex structure is expressed through total abundance and depends on demographic parameters. The resulting one-dimensional equation was applied to describe and estimate demographic characteristics of ungulate populations in the Jewish Autonomous Region. The delay equation provides a good fit to the observed dynamics of ungulate populations, capturing long-term trends in abundance. Point estimates of parameters fall within biologically meaningful ranges and produce population dynamics consistent with field observations. For moose, roe deer, and musk deer, the model suggests predominantly stable dynamics, while annual fluctuations are primarily driven by external factors and represent deviations from equilibrium. Overall, these estimates enable the analysis of structured population dynamics alongside short-term forecasting based on total abundance data.
-
Numerical modeling of flows with flow swirling
Computer Research and Modeling, 2013, v. 5, no. 4, pp. 635-648Views (last year): 4. Citations: 2 (RSCI).This paper is devoted to investigation of the swirl flows. Such flows are widely used in various industrial processes. Swirl flows can be accompanied by time-dependent effects, for example, precession of the vortex core. In turn, the large-scale fluctuations due to the precession of the vortex can cause damage of structures and reduce of equipment reliability. Thus, for engineering calculations approaches that sufficiently well described such flows are required. This paper presents the technique of swirl flows calculation, tested for CFD packages Fluent and SigmaFlow. A numerical simulation of several swirl flow test problems was carried out. Obtained results are compared with each other and with the experimental data.
-
Simulation of corruption in hierarchical systems
Computer Research and Modeling, 2014, v. 6, no. 2, pp. 321-329Views (last year): 8. Citations: 11 (RSCI).Simulation model of corruption in hierarchical systems which takes into account individual strategies of elements and collective behavior of large groups is proposed. Evolution of various characteristics like level of corruption or ratio of corrupted elements and their dependence on external parameters are discussed. The effectiveness of various anticorruptional strategies is examined by means of numeric analysis.
-
Application of simplified implicit Euler method for electrophysiological models
Computer Research and Modeling, 2020, v. 12, no. 4, pp. 845-864A simplified implicit Euler method was analyzed as an alternative to the explicit Euler method, which is a commonly used method in numerical modeling in electrophysiology. The majority of electrophysiological models are quite stiff, since the dynamics they describe includes a wide spectrum of time scales: a fast depolarization, that lasts milliseconds, precedes a considerably slow repolarization, with both being the fractions of the action potential observed in excitable cells. In this work we estimate stiffness by a formula that does not require calculation of eigenvalues of the Jacobian matrix of the studied ODEs. The efficiency of the numerical methods was compared on the case of typical representatives of detailed and conceptual type models of excitable cells: Hodgkin–Huxley model of a neuron and Aliev–Panfilov model of a cardiomyocyte. The comparison of the efficiency of the numerical methods was carried out via norms that were widely used in biomedical applications. The stiffness ratio’s impact on the speedup of simplified implicit method was studied: a real gain in speed was obtained for the Hodgkin–Huxley model. The benefits of the usage of simple and high-order methods for electrophysiological models are discussed along with the discussion of one method’s stability issues. The reasons for using simplified instead of high-order methods during practical simulations were discussed in the corresponding section. We calculated higher order derivatives of the solutions of Hodgkin-Huxley model with various stiffness ratios; their maximum absolute values appeared to be quite large. A numerical method’s approximation constant’s formula contains the latter and hence ruins the effect of the other term (a small factor which depends on the order of approximation). This leads to the large value of global error. We committed a qualitative stability analysis of the explicit Euler method and were able to estimate the model’s parameters influence on the border of the region of absolute stability. The latter is used when setting the value of the timestep for simulations a priori.
-
Simulation of the gas condensate reservoir depletion
Computer Research and Modeling, 2020, v. 12, no. 5, pp. 1081-1095One of problems in developing the gas condensate fields lies on the fact that the condensed hydrocarbons in the gas-bearing layer can get stuck in the pores of the formation and hence cannot be extracted. In this regard, research is underway to increase the recoverability of hydrocarbons in such fields. This research includes a wide range of studies on mathematical simulations of the passage of gas condensate mixtures through a porous medium under various conditions.
In the present work, within the classical approach based on the Darcy law and the law of continuity of flows, we formulate an initial-boundary value problem for a system of nonlinear differential equations that describes a depletion of a multicomponent gas-condensate mixture in porous reservoir. A computational scheme is developed on the basis of the finite-difference approximation and the fourth order Runge .Kutta method. The scheme can be used for simulations both in the spatially one-dimensional case, corresponding to the conditions of the laboratory experiment, and in the two-dimensional case, when it comes to modeling a flat gas-bearing formation with circular symmetry.
The computer implementation is based on the combination of C++ and Maple tools, using the MPI parallel programming technique to speed up the calculations. The calculations were performed on the HybriLIT cluster of the Multifunctional Information and Computing Complex of the Laboratory of Information Technologies of the Joint Institute for Nuclear Research.
Numerical results are compared with the experimental data on the pressure dependence of output of a ninecomponent hydrocarbon mixture obtained at a laboratory facility (VNIIGAZ, Ukhta). The calculations were performed for two types of porous filler in the laboratory model of the formation: terrigenous filler at 25 .„R and carbonate one at 60 .„R. It is shown that the approach developed ensures an agreement of the numerical results with experimental data. By fitting of numerical results to experimental data on the depletion of the laboratory reservoir, we obtained the values of the parameters that determine the inter-phase transition coefficient for the simulated system. Using the same parameters, a computer simulation of the depletion of a thin gas-bearing layer in the circular symmetry approximation was carried out.
-
Numerical simulation of corium cooling driven by natural convection in case of in-vessel retention and time-dependent heat generation
Computer Research and Modeling, 2021, v. 13, no. 4, pp. 807-822Represented study considers numerical simulation of corium cooling driven by natural convection within a horizontal hemicylindrical cavity, boundaries of which are assumed isothermal. Corium is a melt of ceramic fuel of a nuclear reactor and oxides of construction materials.
Corium cooling is a process occurring during severe accident associated with core melt. According to invessel retention conception, the accident may be restrained and localized, if the corium is contained within the vessel, only if it is cooled externally. This conception has a clear advantage over the melt trap, it can be implemented at already operating nuclear power plants. Thereby proper numerical analysis of the corium cooling has become such a relevant area of studies.
In the research, we assume the corium is contained within a horizontal semitube. The corium initially has temperature of the walls. In spite of reactor shutdown, the corium still generates heat owing to radioactive decays, and the amount of heat released decreases with time accordingly to Way–Wigner formula. The system of equations in Boussinesq approximation including momentum equation, continuity equation and energy equation, describes the natural convection within the cavity. Convective flows are taken to be laminar and two-dimensional.
The boundary-value problem of mathematical physics is formulated using the non-dimensional nonprimitive variables «stream function – vorticity». The obtained differential equations are solved numerically using the finite difference method and locally one-dimensional Samarskii scheme for the equations of parabolic type.
As a result of the present research, we have obtained the time behavior of mean Nusselt number at top and bottom walls for Rayleigh number ranged from 103 to 106. These mentioned dependences have been analyzed for various dimensionless operation periods before the accident. Investigations have been performed using streamlines and isotherms as well as time dependences for convective flow and heat transfer rates.
-
Soliton and half-soliton interaction of solitary waves in excitable media with non-linear cross-diffusion
Computer Research and Modeling, 2009, v. 1, no. 1, pp. 109-115Views (last year): 3.We have studied properties of non-linear waves in a mathematical model of a predator – prey system with taxis. We demonstrate that, for systems with negative and positive taxis there typically exists a large region in the parameter space, where the waves demonstrate quasi-soliton interaction; colliding waves can penetrate through each other, and waves can also reflect from impermeable boundaries. In this paper, we use numerical simulations to demonstrate also a new wave phenomenon — a half-soliton interaction of waves, when of two colliding waves, one annihilates and the other continues to propagate. We show that this effect depends on the «ages» or, equivalently, «widths» of the colliding waves.
-
Boundary conditions for lattice Boltzmann equations in applications to hemodynamics
Computer Research and Modeling, 2020, v. 12, no. 4, pp. 865-882We consider a one-dimensional three velocity kinetic lattice Boltzmann model, which represents a secondorder difference scheme for hydrodynamic equations. In the framework of kinetic theory this system describes the propagation and interaction of three types of particles. It has been shown previously that the lattice Boltzmann model with external virtual force is equivalent at the hydrodynamic limit to the one-dimensional hemodynamic equations for elastic vessels, this equivalence can be achieved with use of the Chapman – Enskog expansion. The external force in the model is responsible for the ability to adjust the functional dependence between the lumen area of the vessel and the pressure applied to the wall of the vessel under consideration. Thus, the form of the external force allows to model various elastic properties of the vessels. In the present paper the physiological boundary conditions are considered at the inlets and outlets of the arterial network in terms of the lattice Boltzmann variables. We consider the following boundary conditions: for pressure and blood flow at the inlet of the vascular network, boundary conditions for pressure and blood flow for the vessel bifurcations, wave reflection conditions (correspond to complete occlusion of the vessel) and wave absorption at the ends of the vessels (these conditions correspond to the passage of the wave without distortion), as well as RCR-type conditions, which are similar to electrical circuits and consist of two resistors (corresponding to the impedance of the vessel, at the end of which the boundary conditions are set and the friction forces in microcirculatory bed) and one capacitor (describing the elastic properties of arterioles). The numerical simulations were performed: the propagation of blood in a network of three vessels was considered, the boundary conditions for the blood flow were set at the entrance of the network, RCR boundary conditions were stated at the ends of the network. The solutions to lattice Boltzmann model are compared with the benchmark solutions (based on numerical calculations for second-order McCormack difference scheme without viscous terms), it is shown that the both approaches give very similar results.
-
Method for prediction of aerodynamic characteristics of helicopter rotors based on edge-based schemes in code NOISEtte
Computer Research and Modeling, 2020, v. 12, no. 5, pp. 1097-1122The paper gives a detailed description of the developed methods for simulating the turbulent flow around a helicopter rotor and calculating its aerodynamic characteristics. The system of Reynolds-averaged Navier – Stokes equations for a viscous compressible gas closed by the Spalart –Allmaras turbulence model is used as the basic mathematical model. The model is formulated in a non-inertial rotating coordinate system associated with a rotor. To set the boundary conditions on the surface of the rotor, wall functions are used.
The numerical solution of the resulting system of differential equations is carried out on mixed-element unstructured grids including prismatic layers near the surface of a streamlined body.The numerical method is based on the original vertex-centered finite-volume EBR schemes. A feature of these schemes is their higher accuracy which is achieved through the use of edge-based reconstruction of variables on extended quasi-onedimensional stencils, and a moderate computational cost which allows for serial computations. The methods of Roe and Lax – Friedrichs are used as approximate Riemann solvers. The Roe method is corrected in the case of low Mach flows. When dealing with discontinuities or solutions with large gradients, a quasi-one-dimensional WENO scheme or local switching to a quasi-one-dimensional TVD-type reconstruction is used. The time integration is carried out according to the implicit three-layer second-order scheme with Newton linearization of the system of difference equations. To solve the system of linear equations, the stabilized conjugate gradient method is used.
The numerical methods are implemented as a part of the in-house code NOISEtte according to the two-level MPI–OpenMP parallel model, which allows high-performance computations on meshes consisting of hundreds of millions of nodes, while involving hundreds of thousands of CPU cores of modern supercomputers.
Based on the results of numerical simulation, the aerodynamic characteristics of the helicopter rotor are calculated, namely, trust, torque and their dimensionless coefficients.
Validation of the developed technique is carried out by simulating the turbulent flow around the Caradonna – Tung two-blade rotor and the KNRTU-KAI four-blade model rotor in hover mode mode, tail rotor in duct, and rigid main rotor in oblique flow. The numerical results are compared with the available experimental data.
Indexed in Scopus
Full-text version of the journal is also available on the web site of the scientific electronic library eLIBRARY.RU
The journal is included in the Russian Science Citation Index
The journal is included in the RSCI
International Interdisciplinary Conference "Mathematics. Computing. Education"




