Numerical Modeling of Shockwave Treatment of Knee Joint

Arthritis is a degenerative disease that primarily affects the cartilage and meniscus of the knee joint. External acoustic stimulation is used to treat this disease. This article presents a numerical model of the knee joint aimed at the computer-aided study of the regenerative effects of shockwave treatment. The presented model was verified and validated. A numerical analysis of the conditions for the regeneration of the tissues of the knee joint under shockwave action was conducted. The results allow us to conclude that to obtain the conditions required for the regeneration of cartilage tissues and meniscus (compressive stresses above the threshold value of 0.15 MPa to start the process of chondrogenesis; distortional strains above the threshold value of 0.05% characterized by the beginning of the differentiation of the tissues in large volumes; fluid pressure corresponding to the optimal level of 68 kPa to transfer tissue cells in large volumes), the energy flux density of therapeutic shockwave loading should exceed 0.3 mJ/mm2.


Introduction
Degenerative diseases of the knee joint, such as arthritis and others, significantly reduce a person's quality of life. The average age of the onset of degenerative diseases of this type is 45 years, which leads to disability among a significant part of the working population. However, arthritis can also result from injury [1]. In the end, osteoarthritis has a huge financial toll on the healthcare system [2]. The cartilage plates and meniscus have the function of redefining the loads and helping dissipate the energy obtained as a result of the dynamic effect on the joint. These structural elements of the knee joint are most susceptible to degenerative arthritic changes. That is why the problems of regeneration of cartilaginous plates and meniscus are given a lot of attention in modern medicine. Most acutely, this problem affects the elderly, in whom, due to age, degradation processes begin to appear in the soft tissues and bone tissues of the joints. Former athletes are also highly susceptible to degenerative changes. They have a high level of abrasion of cartilage tissue and soft tissue trauma (meniscus rupture). Surgery is performed in case of catastrophic degenerative changes in the elements of the lower limbs of the skeleton. Pharmacological therapy (based on taking medications, injection, etc.) [3] and non-pharmacological therapy (based on a mechanobiological concept that assumes the division and differentiation of biological cells under the influence of mechanical stimuli) are used to treat mild to moderate degenerative changes in the human musculoskeletal system [4]. Non-pharmacological therapy consists of a set of physical exercises [5] and external physical influence. Shockwave treatment is one of the methods of external exposure therapy [6].
For the treatment of degenerative diseases of the lower extremities, combined therapy is used based on medication (by hormonal and anti-inflammatory drugs) and mechanical treatment (massages, physical exercises, external mechanical stimulation devices). The effect of external mechanical stimulation is based on mechanobiological principles, the essence of which is that a certain level of pressure (mechanical stress) and deformations leads to the growth and differentiation of a certain type of biological tissue. Thus, the regeneration of cartilaginous tissue is facilitated by compressive stresses with an amplitude of 0.15 to 2 MPa (proliferation of chondrocytes) [7]. It was also noted that at stresses lower than 0.003 MPa, chondrogenesis and osteogenesis do not occur, and stresses of the order of compression of 0.7-0.8 MPa are most favorable for the formation of cartilaginous tissue. The optimum for the migration of living cells is the fluid pressure in the pores in the range from 20 kPa to 2 MPa (68 kPa is the most favorable value) [8]. One of the parameters of cartilage tissue differentiation is a value of distortional strain in the range from 0.05 to 1.1% [9,10].
The peculiarities of the physiological state (age, illness, injury) of the patient require an individual selection of the parameters of the mechanical treatment. The energy flux density (EFD) is the main characteristic of the shock wave process for biomedical purposes [11]. The loading is performed through metal or ceramic plates (applicators) of various shapes with a certain level of energy flux density, which is calculated from the parameters of the plates. Currently, active research is underway to determine the optimal loading parameters that contribute to the regeneration of bone and cartilage tissue. Shockwave therapy has been shown to be highly effective in the treatment of knee osteoarthritis [12][13][14]. In the treatment of knee arthritis, applicators are placed in the meniscus and subchondral bone tissue of the tibia [15]. For the treatment of osteoarthritis, shockwave loading with different intensity loading may be effective [16].
Despite the available results, local mechanical changes in the tissues due to acoustic (mechanical) loading of various intensities remain insufficiently understood because of the limitations of experimental studies. The use of computational (in-silico) methods can significantly help in elucidating the mechanical foundations of bone and cartilage tissue regeneration under conditions of low-energy exposure.
The meniscus is the damping area between the femur and the tibia. Due to its special mechanical characteristics, forces are redistributed, and the load on the proximal tibia is reduced [17]. In osteoarthritis, the meniscus is one of the first to change its mechanical characteristics and is also subject to destruction [18]. Therefore, much attention is paid to the study of the conditions for the regeneration of the tissues of the meniscus. To describe the mechanical behavior of meniscus tissues in the nineties of the last century, single-phase viscoelastic models were used [19]. However, such models do not take into account the emerging resistance force from the fluid flow through the pores of the tissues, which, in turn, affects the incorrectness in the description of the mechanical behavior of the tissues under dynamic loading [20][21][22]. In addition, the assumption about the two-parameter condition for the regeneration of biological tissues implies taking into account the liquid. Currently, for studying the mechanical behavior of the knee joint under dynamic action, two-phase poroelastic models are used to describe the material of the meniscus [23]. A cartilage plate performs the antifriction function and promotes energy dissipation in the knee joint [24,25]. It is known that the cartilage tissue of the knee joint has low permeability; this property of the tissue prevents the rapid outflow of the fluid under dynamic loading [26,27]. Therefore, for the numerical study of the conditions of regeneration under dynamic action, two-phase poroelastic models are preferred [28]. Bone tissues may be described by one-phase and two-phase poroelastic models. Poroelastic models are mainly used for the numerical study of the conditions of bone tissue regeneration [29].
Regarding the numerical models devoted to the simulation of the mechanical behavior of the knee joint in the conditions of shockwave therapy, it should be noted that there is no published study that considers the knee as a joint of two bones with articular cartilages and meniscus: there are only the abovementioned publications, which consider different tissues of the joint.
The aim of this work is to develop a numerical model for shock wave exposure on the knee joint as a whole. As a modeling method, the method of movable cellular automata was adopted, which is a representative of computational particle mechanics, and makes it possible to explicitly describe the generation and development of damage in heterogeneous materials. To achieve this goal, a 3D numerical model of the knee joint was verified and validated. Numerical studies were carried out on external shock wave exposure on the knee joint with different energy flux densities.

Method of Movable Cellular Automata
To describe the mechanical behavior of biological tissues, herein we used the model of a poroelastic body, implemented in the method of movable cellular automata (MCA) [30,31], which is an efficient method of computational particle (discrete) mechanics. It has been established that discrete methods have proven themselves to be very promising for modeling contact loading of different materials at both the macroscale and the mesoscale [32,33]. In the MCA method, a solid is considered as an ensemble of discrete elements of finite size (cellular automata) that interact with each other according to certain rules, which, within the particle approach, and due to many body interaction forces, describe the deformation behavior of the material as an isotropic elastoplastic body. The motion of the ensemble of elements is governed by the Newton-Euler equations for their translation and rotation. Within the framework of the method of movable cellular automata, the value of averaged stress tensor in the volume of an automaton is calculated as a superposition of forces that act to the areas of interaction of the automaton with its neighbors [33]. It is assumed that stresses are homogeneously distributed in the automaton volume. Knowing the components of the averaged stress tensor allows adapting different models of plasticity and fractures of classical mechanics of solids to MCA.
The description of the fluid-saturated material in the MCA method is based on the use of effective (implicit) characteristics, such as the volume fraction of interstitial fluid, porosity, permeability, and the ratio of the macroscopic bulk modulus of elasticity to the bulk modulus of the solid skeleton of the material [34]. The fluid filtration in the material is governed by Darcy's law. The mechanical effect of pore fluid on stress and strain of the solid skeleton of the automaton is described using Biot's linear poroelasticity model; therefore, pore fluid pressure affects only the diagonal components of the stress tensor [35].
Details of the method itself and the scheme for numerical solution of the governed equations are given in Appendix A. A three-dimensional version of the method has been implemented in the in-house code MCA3D, which is written in C++ programming language and utilizes Qt library for the pre-and post-processing stages. This code was used in many papers of the authors and their colleagues, in particular, the verification and validation of poroelastic models of tissues of the femur and tibia, based on the MCA method, were carried out in [36][37][38]. A free code for the MCA method solver is also available as a part of the LIGGGHTS package [39]; however, it now provides a reduced functionality for elastic-plastic materials only.

Model of the Knee Joint
The geometric model of the knee joint constructed in this work included the following elements: the epiphyseal part of the femur, the proximal part of the tibia, the cartilage plates, and the meniscus. Each bone consists of a cortical shell and inner cancellous part ( Figure 1).
The knee joint was placed in a box imitating an articular (synovial) capsule, consisting of the interior fibrous tissue and an outer shell (Figure 2a). The applicator was modeled as a thin copper plate of a square shape with a size of 20 × 20 × 1 mm 3 and was located in the meniscus region (black square in Figure 2c). In order to simplify the model and analysis of the results obtained, as well as due to the absence of the need to observe the processes taking place inside the patella, we neglected the presence of the patella. and analysis of the results obtained, as well as due to the absence of the need to observe the processes taking place inside the patella, we neglected the presence of the patella.  The poroelastic properties of the biological tissues of the knee joint adopted in this model are presented in Table 1 and correspond to the data in [38,39]. The fluid in bone tissues is assumed to be equivalent to salt water, with the bulk modulus Kf = 2.4 GPa, and the density ρf = 1000 kg/m 3 [40,41].
In this model, the loading mimicking shockwave exposure was applied using the copper plate as an applicator. To describe the material of the plate, a model of an elastic and analysis of the results obtained, as well as due to the absence of the need to observe the processes taking place inside the patella, we neglected the presence of the patella. The poroelastic properties of the biological tissues of the knee joint adopted in this model are presented in Table 1 and correspond to the data in [38,39]. The fluid in bone tissues is assumed to be equivalent to salt water, with the bulk modulus Kf = 2.4 GPa, and the density ρf = 1000 kg/m 3 [40,41].
In this model, the loading mimicking shockwave exposure was applied using the copper plate as an applicator. To describe the material of the plate, a model of an elastic Standard CAD models available on the Internet (https://www.3dcadbrowser.com/ 3d-model/human-knee-joint, accessed date 1 December 2021) were used as the component parts of the knee joint.
The poroelastic properties of the biological tissues of the knee joint adopted in this model are presented in Table 1 and correspond to the data in [38,39]. The fluid in bone tissues is assumed to be equivalent to salt water, with the bulk modulus K f = 2.4 GPa, and the density ρ f = 1000 kg/m 3 [40,41].
In this model, the loading mimicking shockwave exposure was applied using the copper plate as an applicator. To describe the material of the plate, a model of an elastic body was used with the following parameters: density ρ = 8950 kg/m 3 , bulk modulus K = 115 GPa, shear modulus G = 41.6 GPa [42].
For simulation, we used a computer with Intel i9-10980XE CPU and 64 Gb RAM, running the CentOS 8 operating system. The code MCA3D, implementing the MCA method, utilized parallel computing based on the OpenMP library, so typical computation on 36 threads took about 30-50 h, depending on the number of elements.

Model Verification
The main purpose of verification is to assess the correctness and efficiency of the numerical scheme for solving the governing equations of the method. The key component of numerical model verification is the analysis of the convergence of the obtained results with increasing the resolution of the discrete model (decreasing the discrete element size in the case of computational particle mechanics). The discrete representation of the model is considered optimal when a further increase in its resolution gives no more than a 5% difference with the available resolution [43].
In this work, the analysis for the convergence of a three-dimensional model of the knee joint was carried out by studying the stiffness of the system and the pattern of equivalent strain distribution at different discretization of the considered geometric model ( Figure 2) under its uniaxial compression. Herein, the size (diameter) of discrete elements (automata) in the model sample varied from 0.75 mm to 2.0 mm. The compression of the model samples along the vertical direction was carried out by setting a constant velocity of 0.001 m/s to the upper layer of the particles (Figure 2b).
The results on the convergence of the stiffness of the model knee joint showed that the convergence is nonlinear, and the total scatter between the values for the minimum and maximum size of elements does not exceed 2% ( Figure 3). A very small difference between the values for the size of elements smaller than 1.3 mm indicates a good accuracy of the numerical model for determining its integral parameters. For simulation, we used a computer with Intel i9-10980XE CPU and 64 Gb RAM, running the CentOS 8 operating system. The code MCA3D, implementing the MCA method, utilized parallel computing based on the OpenMP library, so typical computation on 36 threads took about 30-50 h, depending on the number of elements.

Model Verification
The main purpose of verification is to assess the correctness and efficiency of the numerical scheme for solving the governing equations of the method. The key component of numerical model verification is the analysis of the convergence of the obtained results with increasing the resolution of the discrete model (decreasing the discrete element size in the case of computational particle mechanics). The discrete representation of the model is considered optimal when a further increase in its resolution gives no more than a 5% difference with the available resolution [43].
In this work, the analysis for the convergence of a three-dimensional model of the knee joint was carried out by studying the stiffness of the system and the pattern of equivalent strain distribution at different discretization of the considered geometric model ( Figure 2) under its uniaxial compression. Herein, the size (diameter) of discrete elements (automata) in the model sample varied from 0.75 mm to 2.0 mm. The compression of the model samples along the vertical direction was carried out by setting a constant velocity of 0.001 m/s to the upper layer of the particles (Figure 2b).
The results on the convergence of the stiffness of the model knee joint showed that the convergence is nonlinear, and the total scatter between the values for the minimum and maximum size of elements does not exceed 2% (Figure 3). A very small difference between the values for the size of elements smaller than 1.3 mm indicates a good accuracy of the numerical model for determining its integral parameters.   However, the pattern of distortional strain distribution ( Figure 4) indicates sufficient accuracy in determining the zones of strain concentration only for the models with the size of elements of 0.75 and 1.0 mm. Since the calculation times for these models are 13 h and 6 h, respectively, we will take a sample with the diameter of automates equals to 1 mm as the optimal for subsequent calculations.
However, the pattern of distortional strain distribution ( Figure 4) indicates sufficient accuracy in determining the zones of strain concentration only for the models with the size of elements of 0.75 and 1.0 mm. Since the calculation times for these models are 13 h and 6 h, respectively, we will take a sample with the diameter of automates equals to 1 mm as the optimal for subsequent calculations.

Model Validation
The validation of the used poroelastic models of materials for the bone tissues of the knee joint was performed and published in the authors' previous papers [36,37,44].
Validation of the total model of the entire knee joint was carried out by comparing the simulation results with experimental data from [45], and with other numerical simulations from [46]. For this purpose, a compressive load was applied to the upper layer of the femur by its displacement by 0.3 mm for 1 s in accordance with the reference experiment ( Figure 2b).
As mentioned above, the tissues of the cartilaginous plates are subject to the greatest degenerative changes; therefore, the model was validated by comparing the distribution of contact pressure (equivalent stresses in the contact zone) and fluid pressure in these tissues. The distributions of fluid pressure in the pores obtained from our calculations (Figure 5a) are in good qualitative and quantitative agreement with the data presented in the literature [46].
Comparison of the plots for force versus displacement (Figure 5b), which characterizes the rigidity of the system, also showed good agreement with the experimental data presented in [45].

Model Validation
The validation of the used poroelastic models of materials for the bone tissues of the knee joint was performed and published in the authors' previous papers [36,37,44].
Validation of the total model of the entire knee joint was carried out by comparing the simulation results with experimental data from [45], and with other numerical simulations from [46]. For this purpose, a compressive load was applied to the upper layer of the femur by its displacement by 0.3 mm for 1 s in accordance with the reference experiment ( Figure 2b).
As mentioned above, the tissues of the cartilaginous plates are subject to the greatest degenerative changes; therefore, the model was validated by comparing the distribution of contact pressure (equivalent stresses in the contact zone) and fluid pressure in these tissues. The distributions of fluid pressure in the pores obtained from our calculations (Figure 5a) are in good qualitative and quantitative agreement with the data presented in the literature [46]. mm as the optimal for subsequent calculations.

Model Validation
The validation of the used poroelastic models of materials for the bone tissues of the knee joint was performed and published in the authors' previous papers [36,37,44].
Validation of the total model of the entire knee joint was carried out by comparing the simulation results with experimental data from [45], and with other numerical simulations from [46]. For this purpose, a compressive load was applied to the upper layer of the femur by its displacement by 0.3 mm for 1 s in accordance with the reference experiment ( Figure 2b).
As mentioned above, the tissues of the cartilaginous plates are subject to the greatest degenerative changes; therefore, the model was validated by comparing the distribution of contact pressure (equivalent stresses in the contact zone) and fluid pressure in these tissues. The distributions of fluid pressure in the pores obtained from our calculations (Figure 5a) are in good qualitative and quantitative agreement with the data presented in the literature [46].
Comparison of the plots for force versus displacement (Figure 5b), which characterizes the rigidity of the system, also showed good agreement with the experimental data presented in [45]. Comparison of the plots for force versus displacement (Figure 5b), which characterizes the rigidity of the system, also showed good agreement with the experimental data presented in [45].

Simulation of Shockwave Exposure on Knee Joint
The main characteristic of shockwave loading is the value of the energy flux density (PII), which, according to [47], can be expressed through the product of the acoustic wave intensity (I) and normalized time of positive pressure (T p ), as follows: where T p can be defined as the time to reach 90% of maximum positive pressure (Figure 6). At the same time, the intensity (I) is a characteristic of the acoustic impedance of the medium and, in accordance with Equation (1), we obtain an expression for calculating the energy flux density, as follows: where Tp can be defined as the time to reach 90% of maximum positive pressure ( Figure  6). At the same time, the intensity (I) is a characteristic of the acoustic impedance of the medium and, in accordance with Equation (1), we obtain an expression for calculating the energy flux density, as follows: The parameter of "normalized pulse length" was determined from the graphs of pressure versus the time of calculation.
In accordance with the mechanobiological principles, the patterns of the distribution of hydrostatic pressure, the distortional (von Mises) strain, and the fluid pressure in the pores were analyzed. When analyzing the pattern of hydrostatic pressure distribution, it was found that at loading with energy flux density of 0.12 mJ/mm 2 in the area of cartilaginous plates, the conditions for the onset of the processes of osteogenesis and chondrogenesis (hydrostatic pressure higher than 3 kPa) were achieved (Figure 7a). The pattern of distribution of fluid pressure in the pores (Figure 8a  The parameter of "normalized pulse length" was determined from the graphs of pressure versus the time of calculation. In accordance with the mechanobiological principles, the patterns of the distribution of hydrostatic pressure, the distortional (von Mises) strain, and the fluid pressure in the pores were analyzed.
When analyzing the pattern of hydrostatic pressure distribution, it was found that at loading with energy flux density of 0.12 mJ/mm 2 in the area of cartilaginous plates, the conditions for the onset of the processes of osteogenesis and chondrogenesis (hydrostatic pressure higher than 3 kPa) were achieved (Figure 7a Under shockwave loading with an energy flux density of more than 0.3 mJ/mm 2 in the regions of cartilaginous plates, the conditions formed for the regeneration of cartilaginous tissues [7][8][9][10]: the value of compressive stresses exceeds 0.15 MPa (Figure 7b,c), and the values of distortional strain lie in the range from 0.05% to 1% (Figure 8b,c) at values of fluid pressure higher than 20 kPa (Figure 9b,c).
With an increase in the intensity of ultrasonic exposure, the amplitude of the hydrostatic pressure and the area of exposure increase, but the patterns remain practically the same (Figure 7b,c, Figure 8b    With an increase in the intensity of ultrasonic exposure, the amplitude of the hydrostatic pressure and the area of exposure increase, but the patterns remain practically the same (Figures 7b,c-9b,c).
Compressive stresses with maximum amplitude (0.2-1.5 MPa) are concentrated in the region of the cartilaginous plates and meniscus near the loaded surface. With an increase in the amplitude of the energy flux density, the amplitude of the hydrostatic pressure increases, but the pattern remains the same. In the region of maximum stresses, there is also a maximum of values of shear strain up to 0.6% and fluid pressure up to 1 MPa. At the same time, in the bone tissues adjacent to the cartilaginous plate, stresses up to 0.2 MPa are predominantly observed, the shear strain is close to zero, and the fluid pressure in the pores is about 25-30 kPa. This pattern of the distribution of pressure and deformation indicates the creation of conditions for the regeneration of cartilaginous and meniscus tissues, as well as bone tissues near the loaded surface.

Discussion
Shockwave loading of various intensities is used for the regeneration of biological tissues. The acoustic effect is aimed at the fusion of bone tissues (fractures) and the regeneration of tissues subjected to degeneration (cartilage, meniscus). Different types of treatment require a different area of application and power of exposure. It is also necessary to take into account the general degenerative changes in the surrounding bone tissues. There are conflicting data on the effective dosage of shockwave therapy. In addition, ethical restrictions do not allow us to accurately establish the patterns of mechanical stimulation of biological tissues.
In this work, by means of numerical modeling, shockwave loading of various intensities on the knee joint with healthy tissues in the region of the cartilaginous plate of the tibia and meniscus was reproduced. The simulation results indicate the localization of the shock wave action in the area of the loading plate. With an increase in the amplitude of the exposure, the volume in which conditions for tissue regeneration are fulfilled increases. Conditions for chondrogenesis are formed in the cartilage plates, which in turn is consistent with the experimental data on shockwave treatment of the knee joint [15,[48][49][50]. In addition, conditions for chondrogenesis are observed in the tissues of the meniscus, which is also consistent with experimental data [51].
Different authors report conflicting data on the recommended intensity of shockwave loading for the regeneration of knee tissue [12,13]. Earlier, numerical studies on

Discussion
Shockwave loading of various intensities is used for the regeneration of biological tissues. The acoustic effect is aimed at the fusion of bone tissues (fractures) and the regeneration of tissues subjected to degeneration (cartilage, meniscus). Different types of treatment require a different area of application and power of exposure. It is also necessary to take into account the general degenerative changes in the surrounding bone tissues. There are conflicting data on the effective dosage of shockwave therapy. In addition, ethical restrictions do not allow us to accurately establish the patterns of mechanical stimulation of biological tissues.
In this work, by means of numerical modeling, shockwave loading of various intensities on the knee joint with healthy tissues in the region of the cartilaginous plate of the tibia and meniscus was reproduced. The simulation results indicate the localization of the shock wave action in the area of the loading plate. With an increase in the amplitude of the exposure, the volume in which conditions for tissue regeneration are fulfilled increases. Conditions for chondrogenesis are formed in the cartilage plates, which in turn is consistent with the experimental data on shockwave treatment of the knee joint [15,[48][49][50]. In addition, conditions for chondrogenesis are observed in the tissues of the meniscus, which is also consistent with experimental data [51].
Different authors report conflicting data on the recommended intensity of shockwave loading for the regeneration of knee tissue [12,13]. Earlier, numerical studies on acoustic exposure effects have been performed mainly on samples of elementary geometry of biological tissues and biomaterials [52,53]. However, such modeling does not take into account geometric and other physiological features (including the presence of intraosseous fluid) [54] which significantly influence the pressure distribution pattern.
The results of our numerical studies have shown that under shock-wave loading above 0.3 mJ/mm 2 , conditions for the regeneration of the tissues of the knee joint could be achieved. The amplitude of the exposure significantly depends on the volume in which the conditions for tissue regeneration are fulfilled.

Conclusions
The presented numerical model was used for studying shockwave exposure of varying intensity on the knee joint. Analysis of the distribution of hydrostatic pressure, distortional strain, and fluid pressure showed that under shockwave exposure with an energy flux density below 0.3 mJ/mm 2 , conditions for the regeneration of cartilaginous tissues are not formed in the model sample. Under loading with energy flux density higher than 0.3 mJ/mm 2 , the amplitudes of the compressive stress, fluid pressure, and distortional strain corresponding to the conditions of regeneration of cartilaginous tissues are observed in the model samples. Thus, the results obtained are in good agreement with experimental data available from the literature.
In summary, we can conclude that the numerical model of the knee joint developed in this study allows performing a correct simulation of the mechanical behavior of this biological system under shockwave treatment.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.
The method of movable cellular automata (MCA) is a representative of so-called discrete element methods (DEM), which are mainly used for modeling powder and granular materials. The MCA method also possesses some distinguishing features that enable the correct modeling consolidated materials. In this study, the model of poroelastic (two-phase) body was used to describe the mechanical behavior of the biological tissues. Two-phase models are widely used to describe the mechanical behavior of cartilage tissue, but in that case, numerical simulations were only performed based on the computational methods of continuum mechanics (finite elements) [55,56]. Discrete elements were also used to simulate the mechanical response of bone and cartilage tissues but only within the model of the elastic body [57].
In the particle-based method of movable cellular automata [30,31], a simulated body is represented by an ensemble of bonded equiaxial discrete elements of the same size (called movable cellular automata), of which, the position, orientation, and state can change due to interaction with nearest neighbors. Automata interact among each other through their contacts. The initial value of the contact area is determined by the size of the automata and their packing [31,32].
When describing the kinematics and dynamics of an automaton motion, its shape is approximated by an equivalent sphere. This approximation is the most widely used in DEM and allows one to consider the forces of central and tangential interaction of elements as formally independent. This also makes it possible to use the simplified Newton-Euler equations of motion, as follows: where R i , ω i , m i , andĴ i are the position radius-vector, rotation velocity, mass, and inertia moment of i-th automaton, respectively; F pair ij is the pair force of the mechanical interaction of the automata i and j; F Ω i is the volume-dependent force acting on the i-th automaton and is determined by the interaction of its neighbors with other automata. In the second equation of (A1) M ij = q ij n ij × F pair ij + K ij , q ij is the distance from the center of the i-th automaton to the interaction (contact) point of the i-th automaton with the j-th automaton; n ij = (R j −R i )/r ij is the unit vector for orientation of the pair i-j; r ij is the distance between the automata i and j.
For isotropic material, the volume-dependent part of the total force can be written as follows [31,32]: where P j is the pressure in the volume of the neighboring automaton j; S ij is the area of interaction surface of automata i and j; A is defined by the ratio of elastic properties of the material. Let us rewrite the total force acting on automaton i as a sum of normal F n ij and tangential F τ ij components, as follows: where F pair,n ij and F pair,τ ij are the corresponding pair interaction forces depending on the automata overlap h ij and relative tangential displacement l shear ij , respectively. Although the end of Equation (A3) formally corresponds to the interaction force in conventional discrete element models [31], it differs from them essentially due to the many-particle central interaction of the automata.
To compute the components of the average stress tensor in the bulk of automaton i we can use the homogenization procedure described in [31,32]: where α and β denote the global coordinate axes (X, Y, Z); V i is the automaton volume; n ij,α is the α-component of the unit vector n ij ; F ij,β is β-component of the total force acting at the point of contact between automata i and j. The pressure P i (or the mean stress σ i mean ) in the bulk of automaton can be computed from the thus defined stress tensor: We can also compute the other tensor invariants in the bulk of automaton, in particular, the von Mises (equivalent) stress, as follows: The components of the averaged strain tensor ε i αβ are calculated in increments using the specified equation of state of the simulated material and the calculated increments of the averaged stresses. The relation for the force of central interaction of automata is formulated based on the constitutive equation of the material for the diagonal components of the stress tensor, while the force of tangential interaction is formulated on the basis of similar equations for non-diagonal stresses. When implementing the linear elastic model, the expressions for specific values of central and tangential forces of the mechanical response of the automaton i to mechanical action from the neighboring automaton j are written as follows: where the symbol ∆ denotes increment of the corresponding variable for time step ∆t of the numerical scheme of integration of the motion Equation (A1); ∆ε ij and ∆γ ij are the increments of normal and shear strains [31,32] of the automaton i in i-j pair; G i is the shear modulus of the material of the automaton i; K i is the bulk modulus; D i = 1-2G i /K i . Formulas (A1)-(A7) describe the mechanical behavior of a linearly elastic body in the framework of the MCA method. The equations of motion (A1) are usually integrated with the use of the velocity Verlet algorithm, modified by introducing a predictor for estimation of σ i αβ at the given time step [31,32,34]. Due to the necessity of the third Newton's law, the increments of the reaction forces of the automata i and j are calculated based on the solution of the following system of equations: where ∆r ij is the change in the distance between the centers of the automata for a time step ∆t; ∆l sh ij is the value of the relative shear displacement of the interacting automata i and j. The system of Equation (A8) is solved for the increments of strains. This allows calculation of the increments of the specific interaction forces. When solving the system (A8), the increments of mean stress and the values of specific forces in the right-hand sides of relations are taken from the previous time step or are evaluated and further refined within the predictor-corrector scheme.
Automata that model fluid-saturated material are considered as porous and permeable. The pore space of such an automaton can be saturated with fluid. The characteristics of the pore space are taken into account, implicitly, through the specified integral parameters, namely, the porosity φ, the permeability k, and the ratio a = 1 − K/K s of the macroscopic value of bulk modulus K to the bulk modulus of the walls of porous skeleton K s . The mechanical influence of the pore fluid on the stresses and strains in the solid skeleton of an automaton is taken into account on the basis of the linear Biot's model of poroelasticity [28,35]. Within this model, the mechanical response of a "dry" automaton is assumed linearly elastic and is described based on the above-shown relations. The mechanical effect of the pore fluid on the automaton behavior is described in terms of the local pore pressure P pore (fluid pore pressure in the volume of the automaton). In the Biot model, the pore pressure affects only the diagonal components of the stress tensor. Therefore, it is necessary to modify only the relations for the central interaction forces in (A7): ∆F pair,n ij Interstitial fluid is assumed to be linearly compressible. The value of fluid pore pressure in the volume of an automaton is calculated based on of relationships of Biot's poroelasticity model with the use of the current value of pore volume [28,35]. Compressible fluid in the current model is described by the equation of state, as follows: ρ(P pore ) = ρ 0 (1 + (P pore − P 0 )/K fl ), where ρ and P pore are the current fluid density and pressure; ρ 0 and P 0 are the equilibrium values of fluid density and pressure under atmospheric conditions; K fl is the fluid bulk modulus. The pore space of automata is assumed to be interconnected and provides the possibility of redistribution (filtration) of interstitial fluid between the interacting elements. A pore pressure gradient is considered as the "driving force" of filtration. Gravitational effects being neglected, the equation of filtration (fluid mass transfer) in the 'micropore' space can be written as follows: where η is the fluid viscosity and k is the permeability coefficient of the solid skeleton, which can be found as: where d ch is the diameter of the filtration channel. Equations (A9)-(A12) are solved numerically using the Euler scheme (or a higher order numerical scheme) for integration in time on a mesh formed by the centers of the interacting automata (an analog of finite volume method on an ensemble of automata). With regard to the used approximations, there is no mass transfer between the mesh nodes, in which ρ ≤ ρ 0 .