- Open Access
Numerical test concerning bone mass apposition under electrical and mechanical stimulus
Theoretical Biology and Medical Modellingvolume 9, Article number: 14 (2012)
This article proposes a model of bone remodeling that encompasses mechanical and electrical stimuli. The remodeling formulation proposed by Weinans and collaborators was used as the basis of this research, with a literature review allowing a constitutive model evaluating the permittivity of bone tissue to be developed. This allowed the mass distribution that depends on mechanical and electrical stimuli to be obtained. The remaining constants were established through numerical experimentation. The results demonstrate that mass distribution is altered under electrical stimulation, generally resulting in a greater deposition of mass. In addition, the frequency of application of an electric field can affect the distribution of mass; at a lower frequency there is more mass in the domain. These numerical experiments open up discussion concerning the importance of the electric field in the remodeling process and propose the quantification of their effects.
Bones provide mechanical stability to the human body and are a source of minerals for metabolism . Bones have been studied extensively from the mechanical and mineral standpoint, and in terms of functionality [1, 2]. From a mechanical point of view, bones can be adapted to loads on trajectories of stress through mineral apposition, which is due to the action of osteoblasts [1–4]. Furthermore, they reabsorb minerals when the mechanical stimulus is sufficiently low, as it is unnecessary to maintain structure . Reabsorption is directed by osteoclasts. Osteoblasts and osteoclasts are the primary cells involved during bone remodeling that are stimulated by the action of mechanical strain sensors, for example, osteocytes . These three cell types play an important role during the processes of replacement, maintenance and modeling of bones .
Following the work of Meyer during the nineteenth century, Wolff  proposed a theory of trabecular bone architecture, which assumes that trajectories of high mechanical stress form the trabecular bone. In 1987, Frost [6–8] suggested an adaptive mechanism of bone mass as a function of mechanical stress. Consequently, several bone remodeling algorithms have been developed including those proposed by Frost , Pauwels , Kummer , Cowin [11–13] and Hegedus , which predict the formation of bone structure from internal mechanical loads studied in terms of stress and strain.
From mechanical models of bone remodeling, sophisticated studies have been carried out concerning the processes of apposition and reabsorption during bone turnover, and particularly concerning the distribution of mass in the femur [15, 16], hip replacement [17, 18], and dental implants . Generally, these studies have been phenomenological. Therefore, researchers have made significant efforts to include mathematical models and the role of cell biology and biochemistry in the remodeling process, resulting in research that begins at the microscopic level, concerning the effects of basic cellular remodeling units (BMU, Basic multicellular units) during tissue replacement [20, 21]. From the perspective of BMU, important work was initiated at the biochemical and mechanical level concerning the effects of cracks , cell cycles throughout adult life , active molecules within each cell  and the spatial distribution of each BMU . With these important advances in the understanding of bone remodeling, researchers in the field increasingly turned to the study of other biophysical stimuli that can affect this process. Most models have not taken account of the physical-chemical phenomena of tissue mechano-transduction. For this reason, new investigations that allow the piezoelectric and electrokinetic behavior of the bone to be studied were undertaken .
A clinical study demonstrated that a local electromagnetic field accelerates the healing process after bone fracture . Therefore, an article by Demiray and Dost  began new research concerning the effect of the electromagnetic field on interior injury to bone. In another article, Ramtani  presented a mathematical model relating to the benefit of the electric field in the reparair and maintenance of the solid matrix of bone. Furthermore, there are studies concerning the electrical behavior of bone tissue during the production of electric fields, and external electrical flow. Fukada and Yasuda  demonstrated that bone exhibits piezoelectric behavior, i.e. mechanical stress creates electric polarization (the indirect effect) and an external electric field causes strain (the converter effect). In addition, the properties of bones that produce piezoelectric potentials have been determined [29–33]. These data led to the development of mathematical models that include the effect of electromagnetic fields during the repair [34, 35] and remodeling  of bone. For example, Qu and Yu  developed a mathematical model (no spatial dimension) of the remodeling process and healing under the effect of mechanical loads and the use of electric charges. In this model, the higher the voltage applied to a bone after fracture, the lower the percentage of bone damage and micro damage in the few days after the stimulus. Similarly, during osteoporosis an electric field increases bone density over time. Huang et al.  established a hypothesis concerning the biological and biochemical pathways that activate cells, particularly osteocytes, during the imposition of an electric field. Furthermore, Qu and Yu  proposed a mathematical model that included mechanical loads and electromagnetic effects during the process of bone remodeling.
To date, there have been no phenomenological models concerning bone remodeling that have been tested and compared with purely mechanical models. Therefore, this article proposes a new electro-mechanical model relating to bone remodeling. To test its performance, various numerical experiments were carried out and compared with previous mechanical models. The electric model constants were obtained from relevant literature and numerical experimentation. From these assumptions it was concluded that the electric field can affect the distribution of mass, which originates from the remodeling process, under mechanical effects only. Using the remodeling model of Nackenhorst  as a starting point, it was demonstrated that the electric field can increase bone density and accelerate the process of apposition. Therefore, the model proposed herein can be used as a basis for further work concerning electrical effects in the maintenance of bone.
The electro-mechanical model
The electro-mechanical model of bone remodeling that involves mechanical and electrical stimuli can be written, hypothetically, as follows (1):
Where is the well known mechanical stimulus described by Weinans , which depends on tissue density ( ), and the work carried out by mechanical stress ( ) and is the electrical stimulus that depends on density, frequency and the work carried out by the electric field ( ). During this first approach, we consider that the two stimuli are added to determine the bone remodeling process. However, we will develop each of the terms that determine the electro-mechanical model throughout this article.
The mechanical model
Following the remodeling mechanical model described in [4, 39], the density variation over time depends on the mechanical stimulus that exists at every spatial point of the bone, which can be written as in  (2):
Where is the bone density at each point in space ( ), is the energy-strain per unit of volume due to mechanical stress, is a constant and is the strain-energy (per unit of volume) of reference that sets the threshold for which they will perform deposition ( ) or absorption ( ) of the tissue , in the presence of mechanical stress. It should be noted that the energy-strain depends on the density and is given by (3):
Where is the strain, in Voigt notation, of the strain tensor given in (4):
That is a function of the displacements given by (5):
In addition, is the matrix of linear elasticity. The matrix contains the Poisson module, which is usually considered constant, and Young's modulus that depends on the density by expression  (7):
where is a constant and establishes a relationship of power density that has been uncovered through experimentation .
By manipulating equation (7) we can obtain a dimensionless form of the density. This is easier to use with the aim of determining the modulus of elasticity. Multiplying the right side of (7) by produces (8):
Where and are the elasticity modulus and the dimensionless density ratio, respectively. Therefore, the linear elasticity matrix can be expressed as (9):
Where is the matrix of linear elasticity with constant coefficients, which depends on and only, and is given in the case of plane stress by:
Thus, the strain energy per unit of volume (3) can be expressed as (11):
Where is the strain energy at each instant of time, which is calculated with the initial constant of the remodeling  problem only. Replacing these equations in (2), and with some algebraic manipulation, we produce (12):
where we can define . Therefore, we produce the following equation for the density ratio (13):
The momentum equation that establishes the internal stresses of a body is given by  (14):
Where the stress is given by (15):
This article proposes the inclusion of a hypothetical electrical term that can determine, in part, the process of bone remodeling. Therefore, the contribution of this stimulus can be written as (16):
Where is the electrical permittivity of bone tissue and depends on the density, ( ) where frequency ( ), is the electrical energy per unit of volume, is a constant and is the electric energy (per unit volume) of reference.
It should be noted that the term electric energy depends on permittivity, which in turn depends on the density and frequency. This energy term is given by (17):
Where is the electric intensity (electric field). However, , electrical permittivity, depends on density and frequency, and is given by (18):
Where is the permittivity in the free space and is the relative permittivity, which in turn is given by (19):
where is a constant and establishes a relationship of power density that can be proved through experimentation, and will be developed in the following sections. Furthermore, is a function of frequency at which the electric field is applied.
As with the mechanical case, we can manipulate equation (19) so that it is expressed in terms of relative density, multiplying the first term by to produce (20):
where is a constant for the power model of relative permittivity. Therefore, substituting (20) in (17), we obtain (21):
where . Substituting equation (21) in (16) produces (22):
We have chosen as the reference value of electric energy.
However, Gauss's law for an electric field, with no internal loads, is given by (23):
where and is the density of the electric energy. In turn, the electric field can be expressed in terms of a quantity called electric potential or voltage, given by (24):
where is the electric potential.In summary, equation (1) can be written as (25):
Again, using the dimensionless form, we have (26):
where and are mechanical and electrical constants that define the conversion rate of bone remodeling, dependent on the mechanical stress and the electrical potential, respectively.
Solution by the finite element method
Where and are weighting functions, and is the traction on the boundary that serves as a frontier to the domain (see Figure 1).
To calculate the approximate solution by finite element discretization, the displacement field approximates through  (28):
Where is a row vector containing the shape functions used to approximate the displacement  in the nodes. Using the Galerkin method, we approximate the weighting functions in the same way as the displacement field, producing at the elementary level equation  (29):
where is the derivative operator (discrete version of operator of equation (6)) that converts the displacements into strains (see ). This system of equations is completed by applying the Neumann and Dirichlet conditions suitable for solving the elastic problem (see Figure 1).
Similarly, for the electric case we can use the finite element method. For this purpose we use the weak expression of equation (23) in terms of electrical potential to produce (30):
where is the weighting function and is the external normal to the domain of the problem. Again, for calculation of the approximate solution, the scalar field of the electric potential approaches through  (31):
Where is the potential at each node of the element, and are the shape functions for scalar  problems. Weighting functions are chosen in the same manner as shape functions, which is through the method of (Bubnov-) Galerkin standard. Therefore, at the elementary level we have the following equation (32):
where nod is the number of nodes of each element.
Solution to the equation of relative density
This article solves the equation of dimensionless density (or relative) by integrating equation (26) by the method of Euler. For this purpose we define (34):
Therefore, the forward Euler method is defined as  (35):
where is the integration time interval and refers to the evaluation of the variable at a specific time, i.e.: . This method has been used extensively in the prediction of bone density through remodeling . The forward Euler method is of the first order and has the disadvantage of being unstable for large time intervals .
The computational example to be solved is a 0.1x0.1 m square plate with non-uniform load at the top and restrictions of movement vertically at the bottom as presented in Figure 2. This example has been extensively studied in numerical tests concerning bone remodeling algorithms [42–44]. For these values we built a mesh of 80x80 elements (see Figure 3) in the horizontal and vertical directions, respectively, and integrated these with a time step of . The total simulation time was 200. The initial condition was and the limits used in the algorithm of bone remodeling were: .
The constants for the mechanical energy
For the mechanical case we used an initial Young’s modulus , a Poisson's modulus and an initial dimensionless density with values 64 MPa, 0.3 and 1.00, respectively . The dimensionless parameters of the density equation were obtained from Nackenhorst (1997) [39, 44] and were , and .
The mesh was produced using bilinear quadrilateral elements and four points of Gauss integration .
The constants for electric energy
In the electric case, there are numerous articles [45–48] that determine the major electrical properties of bone. In particular, graphs of the electric permittivity in function of frequency and density are provided in , as presented in Figures 4 and 5. Figure 4 presents data concerning bone density as a function of the relative permittivity on the frequency of 50 KHz. The article  used 40 bone tissue samples, and used the measurement of the density as a function of the permittivity. From the data, the mathematical model proposed for permittivity as a function of tissue density is given by a power equation, expressed as (36):
Where and ; this is a function of the density of bone tissue. With the aim of using dimensionless density, we multiplied and divided by to produce (37):
Where and depends on the initial density that considers the computational model.
Figure 5 presents the relative permittivity in terms of the function of the frequency for the density of the femoral head of . With a greater frequency, permittivity decreases. To complement equation (37) (33) we must add a term that allows us to calculate the difference between the relative permittivity at a frequency of 50 kHz and any other frequencies used in the model; this is (38):
where . For its part, the function can be obtained from Figure 5, so we get:
In summary, the equation for the electric permittivity depending on the density and frequency is given by (see Figure 6):
For the remaining constant changes we made variations of and used , as in the mechanical case.
Figure 7 presents various examples that were run for a frequency f = 50 KHz, and with variation in the constant K elect . For these examples, a potential difference at various contour lines was imposed, according to the domain that was established in Figure 2. In the first example, a voltage of 100 was placed on the right side and no voltage at the bottom (Figure 7a). In the second example, a voltage of 100 V was placed on the right side and no voltage on the left (Figure 7b). In the third and fourth examples, the voltage was placed on the top and bottom, respectively (Figures 7c and 7d). In the other contours we placed null Neumann conditions.
Figure 8 presents the results for the relative density value of λ using an approach based on element, with an average at an elementary level of the mechanical stimulus, electrical and density values. In all figures (Figures 8, 9, 10, 11 and 12) black represents the maximum density value (2,175) and white represents the lowest value (0.0125). It should be noted that with Δt = 0.1 the pattern obtained is similar for all values of k elect . It is observed that near the area of imposition of the mechanical stress there is a similar density distribution to a chess board, and far from the loading area there is the formation of three well-defined columns reminiscent of the cortical bone (Figure 8a and 8b). In cases where the constant increases (k elect = 1.4x107 and above) there is the appearance of a high density area in the lower right region, from where a new column starts in the top domain (Figs 8c, d, e and f). In this area, near the lower right structure, there is the distribution and space formation of a chess board. In addition, in the upper part there are empty areas that generate ramifications in each of the columns supporting the load. Note that the imposition of the electric field defines a new topological structure in the domain, as presented on the right-hand side of the simulation results. This new structure is an additional column, which is generated by the potential and a support area of higher density in the lower right corner.
Figure 9 presents the results for various values of K elect . In this case, a voltage of 100 V was imposed on the right side and a null voltage on the left side. Null Neumann conditions were imposed on remaining contours (see Figure 7b). It is noted that there were changes in topology; columns 2 and 3, from left to right in Figures 9c and 9d, are thicker and closer to one another. In addition, upper branches of greater density were created and an additional branch that begins from the last column to the right. In the case of Figure 9e, the formation of a non-defined region of the "chessboard" type on the right side of the domain is presented.
Figures 10 and 11 present results for various values of k elect . In these cases a voltage of 100 V was imposed at the top and bottom, respectively. Null Neumann conditions were imposed on the remaining contours. Note that in these cases the mass distribution increases as k elect increases. In the case of Figure 10e, the density increases in the upper part of the domain, so that the chessboard becomes more continuous in the central part of the domain. In Figures 10e and 11e the columns that are formed in each simulation have higher bandwidths than previous cases.
In the second example there is a variation of the frequency, while the value of k elect remains constant at 7.0x107. With the voltage configuration of Figure 7d, this is with a lower voltage of 100 V. Note in this case that at low frequencies the density and the amount of tissue deposited are greater than at high values. Figure 9a demonstrates that the frequency generates a total deposit of tissue. In Figure 12b the formation of a high density topology and continuing formation is apparent. Figures 12c, 12d and 12e present results comparable to those observed in previous cases. However, in Figure 12c, the columns are wider than in Figures 12d and 12e.
Discussion and conclusions
In this article several numerical examples were developed concerning bone remodeling of the plate during mechanical stress, assuming the imposition of an electric field in the domain. To calculate the mechanical and electrical stimulus of remodeling, and the evolution of density, the elemental approach was utilized. This pioneering article includes the electrical effect, previously designed by Weinans et al. , in a model of bone remodeling.
Comparable with previous articles, the results of the study presented herein demonstrate in the mechanical case formation of the columnar zone (of high density) in the area remote from the load and formation of the trabecular zone (of low density) in the area close to the load . The results are similar to those obtained by Weinans et al. , Fernandez et al.  and Chen et al. . When applying an electric field there is an increase in bone density and an alteration in the topology of the distribution of mass in the domain. In general, there is greater bone mass apposition in the domain. Therefore, the columns developed by the mechanical stress increase in size due to the electric field. In addition, a greater number of columns and localized compact zones may be observed. In the formulation the effect of electrical frequency has been included to allow increased apposition of mass at low frequency to be observed.
Ganong WF, William F: Fisiologia Médica. 2006, Editorial El Manual Moderno, Mexico DF, 373-386.
Cowin SC: Bone mechanics handbook. 2001, CRC Press, Washington
Jacobs CR, Simo JC, Beaupre GS, Carter DR: Adaptive bone remodeling incorporating simultaneous density and anisotropy considerations. J Biomech. 1997, 30: 603-613. 10.1016/S0021-9290(96)00189-3.
Weinans H, Huiskes R, Grootenboer H: The behavior of adaptive bone-remodeling simulation models. J Biomech. 1992, 25: 1425-1441. 10.1016/0021-9290(92)90056-7.
Wolff J: Das Gesetz der Transformation der Knochen. 1892, Hirschwald, Berlin, 1-139.
Frost HM: The laws of bone structure. 1964, CC Thomas, Springfield, IL
Frost HM: Mathematical elements of lamellar bone remodeling. 1964, Thomas, Springfield, IL
Frost HM: Vital biomechanics: proposed general concepts for skeletal adaptations to mechanical usage. Calcif Tissue Int. 1988, 42: 145-156. 10.1007/BF02556327.
Pauwels F: Gesammelte Abhandlungen Zur Funktionellen anatomie Des Bewegungsapparates. 1965, Springer-Verlag, Berlin, 543-549.
Kummer B: Biomechanics of bone: mechanical properties, functional structure, functional adaptation. Biomechanics: its Foundations and Objectives. Edited by: Fung YC, Perrone N, Anliker M. 1972, Prentice-Hall, Englewood Cliffs, 237-272.
Cowin S: Wolffís law of trabecular architecture at remodeling equilibrium. J Biomech Eng. 1986, 108: 83-10.1115/1.3138584.
Cowin S, Hegedus D: Bone remodeling I: theory of adaptive elasticity. J Elast. 1976, 6: 313-326. 10.1007/BF00041724.
Cowin S, Nachlinger RR: Bone remodeling III: uniqueness and stability in adaptive elasticity theory. J Elast. 1978, 8: 285-295. 10.1007/BF00130467.
Hegedus D, Cowin S: Bone remodeling II: small strain adaptive elasticity. J Elast. 1976, 6: 337-352. 10.1007/BF00040896.
Gupta S, New AMR, Taylor M: Bone remodelling inside a cemented resurfaced femoral head. Clin Biomech. 2006, 21: 594-602. 10.1016/j.clinbiomech.2006.01.010.
Stülpner M, Reddy B, Starke G, Spirakis A: A three-dimensional finite analysis of adaptive remodelling in the proximal femur. J Biomech. 1997, 30: 1063-1066. 10.1016/S0021-9290(97)00074-2.
Jonkers I: Relation between subject-specific hip joint loading, stress distribution in the proximal femur and bone mineral density changes after total hip replacement. J Biomech. 2008, 41: 3405-3413. 10.1016/j.jbiomech.2008.09.011.
García J, Doblaré M, Cegoñino J: Bone remodelling simulation: a tool for implant design. Comput Mater Sci. 2002, 25: 100-114. 10.1016/S0927-0256(02)00254-9.
Lian Z: Effect of bone to implant contact percentage on bone remodelling surrounding a dental implant. Int J Oral Maxillofacial Surgery. 2010, 39: 690-698. 10.1016/j.ijom.2010.03.020.
Martin R: Targeted bone remodeling involves BMU steering as well as activation. Bone. 2007, 40: 1574-1580. 10.1016/j.bone.2007.02.023.
Hernandez C, Hazelwood S, Martin R: The relationship between basic multicellular unit activation and origination in cancellous bone. Bone. 1999, 25: 585-587. 10.1016/S8756-3282(99)00201-X.
Taylor D, Tilmans A: Stress intensity variations in bone microcracks during the repair process. J Theor Biol. 2004, 229: 169-177. 10.1016/j.jtbi.2004.03.014.
Hernandez C, Beaupre G, Carter D: A theoretical analysis of the changes in basic multicellular unit activity at menopause. Bone. 2003, 32: 357-363. 10.1016/S8756-3282(03)00037-1.
Peterson MC, Riggs MM: A physiologically based mathematical model of integrated calcium homeostasis and bone remodeling. Bone. 2010, 46: 49-63. 10.1016/j.bone.2009.08.053.
Buenzli PR, Pivonka P, Smith DW: Spatio-temporal structure of cell distribution in cortical Bone Multicellular Units: a mathematical model. Bone. 2010, 48 (4): 918-926.
Ramtani S: Electro-mechanics of bone remodelling. Int J Eng Sci. 2008, 46: 1173-1182. 10.1016/j.ijengsci.2008.06.001.
Demiray H, Dost S: The effect of quadrupole on bone remodeling. Int J Eng Sci. 1996, 34: 257-268. 10.1016/0020-7225(95)00118-2.
Fukada E, Yasuda I: On the piezoelectric effect of bone. J. Phys. Soc. Japan. 1957, 12: 1158-1162. 10.1143/JPSJ.12.1158.
Aschero G, Gizdulich P, Mango F, Romano S: Converse piezoelectric effect detected in fresh cow femur bone* 1. J Biomech. 1996, 29: 1169-1174. 10.1016/0021-9290(96)00011-5.
Beck BR, Qin Y, Rubin C, McLeod K, Otter M: The Relationship of Streaming Potential Magnitude To Strain and Periosteal Modeling 567. Med & Sci Sports & Exercise. 1997, 29: 98-
Gross D, Williams WS: Streaming potential and the electromechanical response of physiologically-moist bone. J Biomech. 1982, 15: 277-295. 10.1016/0021-9290(82)90174-9.
Hung C, Allen F, Pollack S, Brighton C: What is the role of the convective current density in the real-time calcium response of cultured bone cells to fluid flow?. J Biomech. 1996, 29: 1403-1409. 10.1016/0021-9290(96)84535-0.
Johnson MW, Chakkalakal DA, Harper RA, Katz JL: Comparison of the electromechanical effects in wet and dry bone. J Biomech. 1980, 13: 437-442. 10.1016/0021-9290(80)90037-8.
Qu CY, Yu SW: The damage and healing of bone in the disuse state under mechanical and electro-magnetic loadings. Procedia Engineering. 2011, 10: 171-176.
WANG E, ZHAO M: Regulation of tissue repair and regeneration by electric fields. Chinese Journal of Traumatology (English Edition). 2010, 13: 55-61.
Demiray H: Electro-mechanical remodelling of bones. Int J Eng Sci. 1983, 21: 1117-1126. 10.1016/0020-7225(83)90051-4.
Huang CP, Chen XM, Chen ZQ: Osteocyte: the impresario in the electrical stimulation for bone fracture healing. Medical hypotheses. 2008, 70: 287-290. 10.1016/j.mehy.2007.05.044.
Qu C, Qin QH, Kang Y: A hypothetical mechanism of bone remodeling and modeling under electromagnetic loads. Biomaterials. 2006, 27: 4050-4057. 10.1016/j.biomaterials.2006.03.015.
Nackenhorst U: Numerical simulation of stress stimulated bone remodeling. Tech Mech. 1997, 17: 31-40.
Oñate E: Structural analysis with the finite element method. Linear statics. (Springer Verlag, 2009)
Hoffman JD: Numerical methods for engineers and scientists. 1992, Mc Graw-Hill, New York
Beaupre G, Orr T, Carter D: An approach for time dependent bone modeling and remodelingótheoretical development. J Orthop Res. 1990, 8: 651-661. 10.1002/jor.1100080506.
Jacobs CR, Levenston ME, Beaupré GS, Simo JC, Carter DR: Numerical instabilities in bone remodeling simulations: the advantages of a node-based finite element approach. J Biomech. 1995, 28: 449-459. 10.1016/0021-9290(94)00087-K.
Chen G, Pettet G, Pearcy M, McElwain D: Comparison of two numerical approaches for bone remodelling. Medical engineering & physics. 2007, 29: 134-139. 10.1016/j.medengphy.2005.12.008.
Sierpowska J: Prediction of mechanical properties of human trabecular bone by electrical measurements. Physiol Meas. 2005, 26: S119-10.1088/0967-3334/26/2/012.
Sierpowska J: Electrical and dielectric properties of bovine trabecular boneórelationships with mechanical properties and mineral density. Phys Med Biol. 2003, 48: 775-10.1088/0031-9155/48/6/306.
Gabriel C, Gabriel S, Corthout E: The dielectric properties of biological tissues: I. Literature survey. Phys Med Biol. 1996, 41: 2231-10.1088/0031-9155/41/11/001.
Sierpowska J: Interrelationships between electrical properties and microstructure of human trabecular bone. Phys Med Biol. 2006, 51: 5289-10.1088/0031-9155/51/20/014.
Fernández J, García-Aznar J, Martínez R, Viaño J: Numerical analysis of a strain-adaptive bone remodelling problem. Comput Methods Applied Mech Eng. 2010, 199: 1549-1557. 10.1016/j.cma.2010.01.005.
This work was financially supported by Division de Investigación de Bogotá, of Universidad Nacional de Colombia, under the modality: Fortalecimiento a grupos de investigación y creación artística con proyección nacional. ALIANZAS DE GRUPOS. The Project title is “Programa para la obtención y validación de parámetros numéricos utilizados en modelos matemáticos y simulación de sistemas y procesos biológicos mediante el diseño y caracterización de modelos celulares in vitro”, Hermes code 12873.
The authors have no competing interests.
Work concerning production of the manuscript, modelling and numerical simulation was carried out by each of the authors. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.