Smeared Multiscale Finite Element Models for Mass Transport and Electrophysiology Coupled to Muscle Mechanics

Mass transport represents the most fundamental process in living organisms. It includes delivery of nutrients, oxygen, drugs, and other substances from the vascular system to tissue and transport of waste and other products from cells back to vascular and lymphatic network and organs. Furthermore, movement is achieved by mechanical forces generated by muscles in coordination with the nervous system. The signals coming from the brain, which have the character of electrical waves, produce activation within muscle cells. Therefore, from a physics perspective, there exist a number of physical fields within the body, such as velocities of transport, pressures, concentrations of substances, and electrical potential, which is directly coupled to biochemical processes of transforming the chemical into mechanical energy and further internal forces for motion. The overall problems of mass transport and electrophysiology coupled to mechanics can be investigated theoretically by developing appropriate computational models. Due to the enormous complexity of the biological system, it would be almost impossible to establish a detailed computational model for the physical fields related to mass transport, electrophysiology, and coupled fields. To make computational models feasible for applications, we here summarize a concept of smeared physical fields, with coupling among them, and muscle mechanics, which includes dependence on the electrical potential. Accuracy of the smeared computational models, also with coupling to muscle mechanics, is illustrated with simple example, while their applicability is demonstrated on a liver model with tumors present. The last example shows that the introduced methodology is applicable to large biological systems.

Introduction

Convection and diffusion are the fundamental processes regarding transport of molecules or particles (e.g., nanoparticles used in drug delivery) from the cardiovascular system to tissue and interior (cytosol and organelles) (Rushmer, 1976; Hall, 2016). The route of this mass transport goes from blood or lymph to (microenvironment), from extracellular space to intracellular space, and finally to organelles. Along that bidirectional path, molecules/particles pass through the well-known biological barriers: vessel walls and cell and organelle membranes.

The interdisciplinary scientific field where the mass transport is studied, related to drug delivery by applying nanoparticles as the drug carriers to tumor sites, is named oncophysics (Ferrari, 2010; Koay and Ferrari, 2014; Blanco et al., 2015). The fundamental physical principles that govern the mass transport rely on the pressure gradient in case of convection and on concentration gradient in case of diffusion. Of course, the physical process is affected by the physical and chemical characteristics of the transported molecules or particles and the surrounding medium (Reulen et al., 1977; Popel and Pittman, 1998).

Blood as the basic fluid carrying nutrients, oxygen, drugs, and so on is a complex medium, and in computational models, certain approximations have to be adopted. In large vessels, this approximation is made by modeling blood as a homogenous fluid with certain viscosity properties (Shi et al., 2011), where the molecular transport is governed by convective-diffusive laws. This approximation is less valid in small vessels, such as capillaries (diameter size on the order of 10 μm; Mathura et al., 2001) where the presence of cells, mainly red blood cells (RBCs), affects transport (D'Apolito et al., 2015); still, in a simplified analysis, blood may be considered as a homogenous fluid with appropriate viscosity (Sevick and Jain, 1989b; Cinar et al., 1999). Complexity of the capillary network can be seen in various references (e.g., Skinner et al., 1990; Ma and Zhang, 1998), and images are available today on the Internet.

Tissue is a complex medium through which the mass transport occurs, but, regarding the computational modeling, it has to be approximated; the degree of approximation goes from an isotropic continuum to detailed description of extracellular space, cells, and cell interior. Extensive experimental research has been carried out over decades; we cite few recent, as molecule transport within brain tissue (Nicholson, 2001) or (Swabb et al., 1974) where tissue glycosaminoglycan content and drug molecular weight were found to be important parameters determining whether extravascular transport is governed by diffusion or convection. Khaled and Vafai (2003) presented a review of models used for convection and diffusion, including heat transfer, whereas Nugent and Jain (1984) and Gerlowski and Jain (1986) studied diffusivity of dextran molecules in tumor interstitium.

Transport through biological barriers as vessel walls, or cell and organelle membranes, depends on the transport properties of these barriers, as hydraulic or diffusion coefficients, and on the size of the surface that separates the continuum domains. As an illustration, we cite here measurements of the volumes and cross-sectional areas of blood vessels of a dog (Rushmer, 1976), which show that the mass exchange occurs mostly at the capillary–tissue level. The capillary volumetric fraction (capillary density) is one of the important characteristics in the mass exchange (Kojic et al., 2017b), which will represent a basic parameter in our computational models.

As will be seen in the formulation of our computational models, we will rely on the usually adopted governing laws for transport within extracellular space: the Darcy velocity–pressure relationship for convection and Fick's law for diffusion (Kojic et al., 2008). Transport within tumors has additional complexities due to variability of vessel diameters and lengths (Roberts and Palade, 1997) and irregular blood vessel branching (Less et al., 1991), as well as due to geometric resistance (Sevick and Jain, 1989a) (a measure of network irregularities), viscous resistance (Sevick and Jain, 1989b), and RBC mechanical properties (Sevick and Jain, 1991). Jain (1988) summarized basic properties of blood flow within tumors, while data specifically related to transport through capillary walls are given in Jain (1987).

The fundamental process in living organisms represents mass exchange and transport within cells. Cell interior is very complex, composed of different compartments or entities, such as organelles, which are, together with the cell cytoskeleton, immersed into the cell fluid–cytosol (Keener and Sneyd, 2009). Mass transport within cells is a particular branch of biomedical research, and here we cite a few references addressing some aspects of mass transport and exchange. The factors that affect the transport within cells range from biochemical to mechanical to signaling pathways (Rangamani and Iyengar, 2007). Theoretical basis of biochemical processes within different intracellular compartments is formulated (Lauffenburger and Horwitz, 1996; Sun and Pang, 2008; Chu et al., 2013). A computational framework for modeling mass exchange within cells, formulated as a virtual cell, is presented (Schaff et al., 1997; Moraru et al., 2008), which was further used in many applications (Slepchenko et al., 2003).

In electrophysiology, the goal is to determine the electrophysiological properties of all compartments and signal propagation characteristics within the body. For example, in heart electrophysiology, the fundamental advancements were achieved by designing the so-called clamp experiment (Hodkin and Huxley, 1952) in order to determine characteristics of the membrane currents and constitutive relations for the conduction of currents. Following this breakthrough achievement, further important experiments were performed, with modifications and extension of the constitutive relations (Noble, 1962; Baer and Rinzel, 1991; Decker et al., 2009; O'Hara et al., 2011). Based on these experimental results, a large number of computational models related to electrical signal propagation, which includes enormous complexities, have been introduced. For example, model presented (Noble, 1962) includes the role of various transmitting molecules—currents carried by ions through membranes and composite media (Ijiri et al., 2008) include intricate geometry such as in the case of Purkinje network. Since heart is a representative organ with electrical signaling coupled to mechanics, we here refer to a few computational models particularly related to heart modeling and mainly based on the finite element method. A review of computational methods for heart physiology can be found (Clayton and Panfilov, 2008). Among models discussed (Clayton and Panfilov, 2008) is the one related to coupling the electrical field with the mechanical response (Rocha et al., 2011; Lafortune et al., 2012; Dal et al., 2013).

Finally, in this short overview of mass transport and electrophysiology fields of research in science, technology, and medicine, we note that variation of electrical potential of cell membranes triggers other vital processes within living cells. A typical example is that calcium waves within muscle cells are induced by changes of the membrane potential. Calcium concentration within cells is fundamental for muscle contraction, and various mechanical models for muscles based on the calcium concentration changes are developed (Hunter et al., 1998; Kim et al., 2010; Mijailovich et al., 2010; Lafortune et al., 2012; Santiago, 2018). The variation of electrical potential fields also affects transport of charged particles or drug molecules (Schaff et al., 1997; Trapp and Horobin, 2005), and this kind of mass transport is included in our models (Kojic et al., 2019).

Considering the complexity of computational modeling of physical fields, such as pressure distribution, concentration of various molecules, or field of electrical potential within a body, or even within an organ, it is desirable to have a computational methodology feasible for practical applications. In our previous research, we have introduced a smeared concept for mass transport in capillary system and tissue (Kojic et al., 2017a,b, 2018, 2019; Kojic, 2018; Milosevic et al., 2018a) and demonstrated its superiority with respect to traditional modeling methods.

In this study, we summarize our previously formulated smeared methodology for modeling mass transport in biological tissue and electrophysiology problems and extend it to coupling with muscle mechanics. The motivation for the development of the models, which are applicable to real biomedical problems, is briefly outlined in the text above.

The paper is organized as follows. In the next section, we first summarize the governing equations of the gradient-driven processes: fluid flow through porous media, convective–diffusive particulate transport (including ions) and electrostatics; and then introduce a finite element (FE) form for these governing equations. In Section Smeared Model for Field Problems, the smeared methodology in a general form is summarized (Kojic, 2018), followed by specific forms related the problems described in Section A Summary of the Fundamental Equations for Gradient-Driven Physical Processes and FE Formulation. In Section Coupling Electrophysiology and Muscle Mechanics, we give a short description of coupling the smeared model for electrical potential with muscle mechanics. In Section Representative Results, we present several examples to demonstrate the main features of the smeared concept, including its accuracy and applicability; and, in the final section, we give a short summary and concluding remarks.

A Summary of the Fundamental Equations For Gradient-Driven Physical Processes and FE Formulation

In this section, we first summarize the gradient-driven problems related to mass transport in blood vessels and tissue and electrophysiology. Then, we present a finite element formulation for these partial differential equations.

Fundamental Equations for the Gradient Driven Field Problems

Flow Through Porous Media

In case of incompressible fluid flow through a porous rigid medium, the governing relation is represented by Darcy's law

$\begin{array}{ll}{v}_i=-k_{D_{ij}}\frac{\partial p}{\partial x_j},sumonj:j=1,2,3& (1)\end{array}$

where v_i is the Darcy velocity (as fluid flux per unit area of the continuum) in direction x_i, p is fluid pressure, and k_Dij is the Darcy tensor. The mass balance equation is

$\begin{array}{ll}{k}_{D_{ij}}\frac{{\partial }^{2}p}{\partial x_i\partial x_j}+q_V=0& (2)\end{array}$

WHERE q_V is a source term.

Diffusion

The constitutive law for diffusion is known as Fick's law

$\begin{array}{ll}{Q}_i=-D_{ij}\frac{\partial c}{\partial x_j}& (3)\end{array}$

and the mass balance equation is

$\begin{array}{l}-\frac{\partial c}{\partial t}-v_i\frac{\partial c}{\partial x_i}+\frac{\partial }{\partial x_i}\left(D_{ij}\frac{\partial c}{\partial x_j}\right)+q_V=0,\\ sumoni,j:i,j=1,2,3& (4)\end{array}$

Here, c is concentration, Q_i flux, and D_ij is the diffusion tensor. The generality is kept under the assumption that the diffusion tensor can be a function of concentration, that is, it can be D_ij = D_ij(c).

Electrostatics

The constitutive law is

$\begin{array}{ll}{J}_i=-G\frac{\partial V_e}{\partial x_i}& (5)\end{array}$

where G is electric conductivity, and V_e is electrical potential. The continuity equation for the current density can be derived from Maxwell equations in the form

$\begin{array}{ll}\frac{\partial }{\partial t}\left(\frac{\partial D_i}{\partial x_i}\right)=-\frac{\partial J_i}{\partial x_i},sumoni,i=1,2,3& (6)\end{array}$

where the current density components D_i can be related to the potential V_e as

$\begin{array}{ll}{D}_i=-\epsilon \frac{\partial V_e}{\partial x_i}& (7)\end{array}$

where ε is the dielectric constant. Finally, the continuity equation is

$\begin{array}{ll}\epsilon \frac{\partial }{\partial t}\frac{{\partial }^2{V}_e}{\partial x_i\partial x_i}=-G\frac{{\partial }^2{V}_e}{\partial x_i\partial x_i}+q_V^e& (8)\end{array}$

where $q_V^e$ is a volumetric source term (coming from ionic transport; Rangamani and Iyengar, 2007).

1D Conditions

For further presentation, we give the expressions for the 1D conditions. For the fluid flow, the 1D conditions follow from the study of flow within pipes (Smith et al., 2002; Canic and Kim, 2003; Kojic et al., 2014). In case of a rigid pipe, the governing equation reduces to

$\begin{array}{ll}{k}_{pipe}\frac{{\partial }^{2}p}{\partial {\overline{x}}^2}=0& (9)\end{array}$

where $\bar{x}$ is the pipe direction, and k_pipe is the pipe coefficient, which can be derived from the so-called Hagen-Poiseuille law. Additional terms are present in the above equation for the case of deformable pipe (Kojic et al., 2014) but will not be considered in this work.

In case of diffusion, the 1D conditions follow from Equation (4). Hence, we have

$\begin{array}{ll}-\frac{\partial c}{\partial t}-v\frac{\partial c}{\partial \overline{x}}+\frac{\partial }{\partial \overline{x}}\left(D\frac{\partial c}{\partial \overline{x}}\right)+q_V=0& (10)\end{array}$

where $\bar{x}$ is the axis of propagation, and D is diffusion coefficient. In the case of electrical conduction, the governing equation has the form (9) with respect to the electric potential V_e, where instead of k_pipe, we have G_aA, with A being the neural fiber cross-section.

Transport Through Membranes

Continuum domains of a composite media are often separated by membranes or walls in case of blood vessels or neural fibers. For the presentation of the smeared methodology, we here give the fundamental relations for transport through membranes. In case of fluid flow or diffusion, we have

$\begin{array}{ll}{Q}_w^p=k_w\left(p_{in}-p_{out}\right)& (11)\end{array}$

$\begin{array}{ll}{Q}_w^c=D_w\left(c_{in}-c_{out}\right)& (12)\end{array}$

with the flux of fluid $Q_w^p$ and mass due to diffusion $Q_w^c$ oriented outward (from in to out); k_w and D_w are wall hydraulic permeability and wall diffusivity, respectively. In the case of electrical field, the wall electrical flux relies on the so-called cable theory (Winslow, 1992). The outlet electrical flux (current density) I_m can be expressed as

$\begin{array}{ll}{I}_m=G_m\left(V_e^{in}-V_e^{out}\right)+C_m\left(\frac{\partial V_e^{in}}{\partial t}-\frac{\partial V_e^{out}}{\partial t}\right)& (13)\end{array}$

where G_m is membrane conductivity, and C_m is specific membrane (wall) capacitance.

Finite Element Formulation

The above governing equations can be transformed into the FE equations of balance for a single finite element by implementing a standard Galerkin weighting method (Kojic et al., 2008). The incremental-iterative form of balance for a time step Δt and iteration i can be derived in the form

$\begin{array}{l}{\left(\frac{1}{\Delta t}M+K^v+K\right)}^{\left(i-1\right)}\Delta {\Phi }^{\left(i\right)}=Q^{ext}+Q^V\\ -\frac{1}{\Delta t}{M}^{\left(i-1\right)}\left({\Phi }^{\left(i-1\right)}-{\Phi }^t\right)-{\left(K+K^v\right)}^{\left(i-1\right)}{\Phi }^{\left(i-1\right)}& (14)\end{array}$

where

$\begin{array}{l}{M}_{IJ}={\displaystyle \underset{V}{\int }{c}_m{N}_I{N}_J}dV\\ K_{IJ}^v={\displaystyle \underset{V}{\int }{v}_i{N}_I}{N}_{J,i}dV,sumoni:i=1,2,3\\ K_{IJ}={\displaystyle \underset{V}{\int }{D}_{km}{N}_{I,k}{N}_{J,m}}dV,sumonk,m:k,m=1,2,3\\ Q_I^V={\displaystyle \underset{V}{\int }{N}_I{q}_{V}dV}& (15)\end{array}$

where Φ stands for pressure, concentration, or electrical potential as nodal variables; N_I are interpolation functions, V is element volume; c_m is mass coefficient (= 0 for fluid flow, and = 1 for diffusion); D_km for fluid is the Darcy tensor, while it is Gδ_ij (δ_ij is the Kronecker delta symbol) for electrical field. For the case of Darcy's flow or no convection, the convection matrix K^v is equal to zero. For electrical potential, we have that the “mass” matrix is

$\begin{array}{ll}{M}_{IJ}=\epsilon {\displaystyle \underset{V}{\int }{N}_{I,k}{N}_{J,k}}dV,sumonk,k=1,2,3& (16)\end{array}$

and the convection matrix is equal to zero. Note that for 1D problems, the equations have the same form as the above, with one index k and no summation; and that the element volume is V = AL, where A is cross-sectional area and L is the element length. Note that Equation (14) assumes implicit integration scheme over time, that is, all variables are evaluated at the end of time step and at the current equilibrium iteration. This integration scheme is unconditionally stable and provides the best accuracy (Kojic and Bathe, 2005).

For modeling transport through the membranes (walls), we have introduced connectivity 2-node elements for nodes at membranes (Figure 1A) by using double nodes at the same space position at the membrane, with one node belonging to one side and the other to the other side of the boundary between two domains. The balance equation of the form (14) can be applied. The “mass” and transport matrices M and K can be written as

$\begin{array}{l}{M}_{11}=M_{22}=\frac{1}{3}{c}_{mm}{A}_m{h}_m,M_{12}=M_{21}=\frac{1}{6}{c}_{mm}{A}_m{h}_m\\ K_{11}=K_{22}=-K_{12}=-K_{21}=D_w{A}_m& (17)\end{array}$

where A_m is the area of the surface belonging to the node, h_m is the membrane thickness; in case of diffusion, D_w is the diffusion coefficient, while instead of D_w, we have k_w and G_m for fluid flow and for electrical conduction, respectively; c_mm = 0 for fluid flow and c_mm = 1 for diffusion. For the case of electrical conduction, the non-zero terms of the “mass” matrix are

$\begin{array}{ll}{M}_{11}=M_{22}=A_m{C}_m& (18)\end{array}$

FIGURE 1

Figure 1. Schematic of detailed model and smeared model. (A) Detailed model of tissue as composite medium with continuum subdomains and capillaries/fibers, 2D representation, with continuum, 1D and connectivity elements. (B) Smeared finite element (FE) representation of the detailed model. (C) Composite smeared finite element (CSFE) with subdomains and connectivity element at a FE node J. (Used with permission from the Journal of the Serbian Society for Computational Mechanics from Kojic, 2018).

Finally, in this section, we summarize computational steps related to the incremental-iterative balance Equation (14). The steps are as follows:

1. Loop over time steps

1.1 Loop over equilibrium iterations

1.1.1 Loop over finite elements, evaluate, and store into global system element matrices and nodal vectors

1.1.2 Solve for the nodal variables

1.1.3 Check for the convergence criteria (unbalanced vector on the right-hand side are small and/or increments of nodal variables are small).

If the convergence criteria are not satisfied, update the solution and go to next equilibrium iteration

1.2.1 If convergence criteria are satisfied, store solution and go to the next time step 1.

END OF SOLUTION.

Smeared Model for Field Problems

To introduce the smeared methodology (Kojic et al., 2014), we first consider a detailed model of a composite medium. A schematic of a medium composed of continuum domains-compartments together with a network of fiber-like 1D domain is shown in Figure 1A. The continuum domains include extracellular space, cells, and organelles. Capillaries and lymph vessels, and neural fibers, are represented by 1D elements within extracellular space, while cells can contain organelles—hence, the continuum domains have a hierarchical character. Each domain has its own FE mesh of continuum elements, while 1D domains have their own 1D finite elements with the coordinate axes along the elements ($\bar{x}$ axis depicted at one of these fibers).

Connectivity elements are shown as A, B; C, D; E, F; and enlarged at the top of the Figure 1A (the nodes C, D, also denoted as 1, 2). The two nodes, 1 and 2, have nodal values representing the two domains (ϕ_out and ϕ_in in the figure). The characteristics of the connectivity elements assigned to the boundary common nodes are membrane (or wall) transport coefficient, cross-section area A_m, and the membrane (wall) thickness h_m.

It can be seen from the detailed model description that significant effort is required to generate the model. The model generation would be an impractical or even impossible task in case of a complex medium such as tissue. In case of employing continuum FEs for membranes instead of connectivity elements, this task would be more demanding.

To overcome complexity issues, we further introduce a smeared model by formulating a continuum composite finite element (CSFE). The CSFE includes all constituents (continuum and 1D) in a way that the true physical fields, corresponding to a detailed model, are represented in a kind of average, so-called “smeared” sense, providing adequate accuracy. The smeared model for the same detailed model of Figure 1A, and with continuum elements only, is schematically shown in Figure 1B.

A few conceptual steps need to be considered to formulate the CSFE element. First, the 1D balance equations are necessary to be transformed into the corresponding continuum format. The derivations of the continuum transport tensor, for a general physical field, are given (Kojic, 2018) and can be expressed as

$\begin{array}{ll}{D}_{ij}=\frac{1}{A_{tot}}{\displaystyle \sum _K{D}_K{A}_K{\ell }_{Ki}}{\ell }_{Kj}& (19)\end{array}$

where A_tot is the total area of 1D domains in a reference volume; A_K is cross-sectional area of element K; D_K and ℓ_Ki are transport coefficients along the 1D elements and directional coefficients, respectively.

Each domain in the CSFE formulation has its own field within the corresponding volume of the CSFE. Therefore, as shown in Figure 1C, the FE node of the CSFE has a number of degrees of freedom ϕ^K equal to the number of domains N_d. The domain volume V_K can be expressed as

$\begin{array}{ll}{V}_K=r_V^{K}V,\text{ and d}{V}_K=r_V^{K}dV& (20)\end{array}$

where $r_V^K$ is the volumetric fraction, and V is the total element volume.

Corresponding domains are coupled using connectivity elements placed at each node of the CSFE and in accordance to previously described connectivity elements in the detailed model. The cross-sectional area A_JK of a connectivity element at node J for the domain K can be expressed in the form

$\begin{array}{ll}{A}_{JK}={\left(r_{AV}^K{V}_K\right)}_J={\left(r_{AV}^K{r}_V^{K}V\right)}_J& (21)\end{array}$

where $r_{AV}^K$ is the area coefficient expressed as $r_{AV}^K=A_K/V_K$, and V_J is the volume of the total space of the continuum belonging to the node J. This volume can be evaluated as

$\begin{array}{ll}{V}_J={\displaystyle \sum _{el}{\displaystyle \underset{V_{el}}{\int }{N}_{J}d{V}_{el}}}& (22)\end{array}$

where summation el goes over all finite elements with the common node J. Note that for convenience of modeling of any non-homogenous property of the material system, all of the surfaces, volumes, and volumetric and area ratios are assigned to FE nodes.

The finite element balance equations for continuum and connectivity elements are of the same form as in detailed model (Equation 14) with the matrices

$\begin{array}{l}{M}_{IJ}={\displaystyle \underset{V}{\int }{c}_m{N}_I{N}_J{r}_V^K}dV\\ K_{IJ}={\displaystyle \underset{V}{\int }{D}_{ij}{N}_{I,i}{N}_{J,j}{r}_V^K}dV,sumoni,j:i,j=1,2,3\\ Q_I^V={\displaystyle \underset{V}{\int }{N}_I{q}_V{r}_V^{K}dV}& (23)\end{array}$

where the material parameters are as in Equation (15). In case of electrical field, the matrix M in Equation (16) is now

$\begin{array}{ll}{M}_{IJ}=\epsilon {\displaystyle \underset{V}{\int }{N}_{I,k}{N}_{J,k}}{r}_V^{K}dV,sumonk,k=1,2,3& (24)\end{array}$

and the source nodal vector due to ionic transport of a molecule m is (Kojic et al., 2019)

$\begin{array}{ll}{Q}_I^{mE}=\frac{DF{z}^m}{RT}{\displaystyle \underset{V}{\int }{N}_I\frac{\partial }{\partial x_i}\left(c^m\frac{\partial V_e}{\partial x_i}\right)}{r}_V^{K}dV& (25)\end{array}$

where D is diffusion coefficient, z^m is the molecule valence, F is the Faraday constant, R is the gas constant, T is absolute temperature, and c^m is concentration. The surface areas entering into matrices of the connectivity elements are as given in Equation (21).

The above concept has been implemented for diffusion and fluid transport through complex systems consisting of capillary network and tissue, with inclusion of the hydrophobicity effects in the connectivity elements (Kojic et al., 2014, 2017a, 2018), and with improvements of accuracy of the smeared methodology achieved by the correction function introduced in Milosevic et al. (2018a). The robustness and applicability of the smeared model are demonstrated (Milosevic et al., 2018b), where the CSFE is extended for modeling drug release from a complex mesh of drug-loaded nanofibers.

Coupling Electrophysiology and Muscle Mechanics

Muscles (here assumed skeletal muscles) in the body are activated by electrical signals transmitted from the central nervous system to muscle cells. The signals trigger muscle activation since they produce a change in the cell membrane potentials, which further leads to flow through membrane of ions vital for cell functioning, such as potassium, sodium, calcium, and others (Hodkin and Huxley, 1952; Noble, 1962; O'Hara et al., 2011). The ion flow is bidirectional through various biological mechanisms. There are a number of mathematical models that connect the membrane potential change with activation of muscles. For example, for cardiac muscle, the mathematical expressions for generation of the so-called active stress along the muscle fiber, which produces the muscle contractile force, the membrane potential is used directly (Dal et al., 2013) or through the concentration of calcium Ca²⁺ within the muscle cell (Hunter et al., 1998; Lafortune et al., 2012; Berberoglu et al., 2014). The calcium Ca²⁺ is the crucial molecule that catalyzes the biochemical cycle of conformational change of muscle fiber molecules and therefore transformation of chemical into mechanical energy. Hence, in modeling muscle mechanical action, it is necessary for these models to have the calcium concentration change within muscle cells over time. We will use a widely accepted relation (Hunter et al., 1998)

$\begin{array}{ll}{\sigma }_{act}=\frac{{\left[C{a}^{2+}\right]}^n}{{\left[C{a}^{2+}\right]}^n+C_{50}^n}{\sigma }_{\max }\left[1+\eta \left(\lambda -1\right)\right]& (26)\end{array}$

where σ_act is the active stress along the fiber, Ca²⁺is calcium concentration, σ_maxis maximum isometric stress, C₅₀ is concentration for 50% availability of actin sites for the crossbridge binding, n is related to the rate of this availability to concentration, η is parameter that is governing the rate of muscle fiber deformation, and λ is the fiber stretch.

Computational models of skeletal muscles and their implementation to the physiological conditions have been the subject of intensive research within the computational community (e.g., Kojic et al., 1998; Stojanovic et al., 2007; Mijailovich et al., 2010, and references therein).

We use velocity formulation, that is, the nodal variables are velocities—convenient to couple solid and fluid mechanics (Isailovic et al., 2014; Kojić et al., 2015), while stresses in solids are calculated from strains or stretches. The balance equation of a finite element can be written in the form (Kojic et al., 2008)

$\begin{array}{ll}\left(\frac{1}{\Delta t}M+\Delta tK\right)\Delta V^{\left(i\right)}=F^{ext}-F^{\mathrm{int}\left(i-1\right)}-\frac{1}{\Delta t}M\left(V^{\left(i-1\right)}-V^t\right)& (27)\end{array}$

where the mass and stiffness matrices M and K have a standard form (Kojic et al., 2008), and V and V^t are nodal velocities at the current (or previous) iteration and at start of time step, respectively; F^ext and F^int are external and internal nodal forces; nodal variables are one-dimensional arrays. Specific to muscle deformation is that, besides the material stress dependent on the state of deformation, there exists the active stressσ_act entering into the internal force vector as noted above.

Representative Results

We have selected two numerical examples to illustrate the main characteristics of our smeared modeling methodology. The first example shows applicability of smeared models to large biological systems; here, a liver with two tumors is modeled. In the second example, accuracy of solutions is investigated on a small sample of the heart wall tissue, including coupling between electrophysiology and mechanics. The models are built into our FE code PAK (Kojic et al., 2010) with the corresponding CAD interface developed at R&D Center for Bioengineering BIOIRC.

We prepared animations for both examples, which are given in Supplementary Material.

Liver Model With Tumors—Application of the CSFE

Liver model is based on geometry used (Kojic et al., 2017b), with two additional spheres modeled to mimic two tumor domains. Micro computed tomography (micro-CT) was used by the Preclinical Imaging Core at the Houston Methodist Research Institute to scan the vascular structure of a mouse liver. The geometry of the liver and blood vessel network is generated from micro-CT scan using procedure presented (Zagorchev et al., 2010). The geometry of the model is shown in Figure 2.

FIGURE 2

Figure 2. Liver model: geometry, tumor domains, and pressures within large vessels.

The FE model consists of 39,832 3D CSFEs, 7,736 1D pipe FEs for larger vessels, 726 connectivity elements for connecting large vessels with continuum nodes of smeared FEs. Two separate tumors within liver are with a total of 316 3D FE elements. The total number of nodes is 54,590. Data for this example are the same as in Kojic et al. (2017b).

Prescribed conditions in larger vessels (at input/output nodes of 1D pipe elements mesh) are: Inlet Pressure is 3,999.7 Pa (30 mmHg), Outlet Pressure is 1,333.2 Pa (10 mmHg), Inlet Concentration is bolus type c(t), presented in Figure 6, and outlet concentration is set to be 0.

Characteristics of fluid/diffusion flow through blood vessels (large vessels and capillaries) are: viscosity is 10⁻³ Pa s, and diffusion coefficient is 1,000 mm²/s.

Characteristics of blood vessel walls are: hydraulic permeability coefficient is 10⁻¹² mm/(Pa s), diffusion coefficient is 0.1 mm²/s, and thickness is 10% of the vessel diameter.

Characteristics of tissue are: diffusion coefficient is 0.1 mm²/s, and Darcy coefficient is 10⁻¹² mm²/(Pa s).

Parameters of smeared model are: average capillary diameter is 0.025 mm, capillary wall thickness is 0.0025 mm, volume fraction is 10%, diffusion coefficient of capillary wall is 10⁻⁶ mm²/s, and hydraulic permeability of capillary wall is 10⁻¹² mm/(Pa s). The adopted values of material and geometric parameters are according to literature data, for example, Rushmer (1976), Gerlowski and Jain (1986), Jain (1988), Keener and Sneyd (2009).

Time steps used in simulations are: 40 × 2.5 s.

Additionally, there are two tumors in the model with the following characteristics:

As can be seen from these data, we assumed 10 times smaller diffusion coefficient within tumor 2 than in tumor 1 to show the difference in concentration between these tumors.

With these material data and boundary conditions, we have solved for pressures and concentrations within the liver and two tumors using our smeared methodology. Some of the results are summarized below.

Pressure fields for two views of the model are shown in Figure 3 for the outer surface of 3D smeared elements, cross-section, and a dotted representation of large vessels and continuum. Tumor surfaces are indicated by dashed lines. As can be seen from this figure, there is evident reduction in pressures, starting from large vessels to capillaries and further to tumors and healthy tissue.

FIGURE 3

Figure 3. Pressure fields for two views: (A) full mesh; (B) clipped mesh; (C) dotted representation of results in tissue, and with full mesh in tumors.

Vectors of velocities, field of velocity on the outer surface with a clipping plane, and vectors of velocities within tumors, for two different views on the model, are shown in Figure 4. For one of the tumors, T1, there is a positive pressure difference with respect to the surrounding tissue, which induces velocity vectors pointing out of tumor, that is, fluid is flowing out of tumor (Figure 4C–T1).

FIGURE 4

Figure 4. Velocities within liver model for two views: (A) vectors of velocities; (B) outer surface with clipped mesh; (C) vectors of velocity within tumors (enlarged images).

The concentration field within large vessels, liver tissue, and tumors is shown in Figure 5 for three time points. The largest values of the concentrations are noticeable in blood vessels, following a decrease going to capillaries and tissue. Also, concentration within tumor T2 is smaller compared to that of tumor T1 due to reduced diffusion coefficients.

FIGURE 5

Figure 5. Concentration field in liver with tumors (marked with dashed lines), dotted results in tissue domain and with full mesh in tumors, for times t = 10, 20, and 50 s.

Change of the mean concentration within capillaries, tissue of liver, and tumors T1 and T2, is shown in Figure 6. In the first period of the process of drug transport within the liver, concentration within capillaries, tissue, and tumors is increasing according to increase of the entering mass (concentration). The concentration within capillaries reaches maximum with respect to the maximum of the entering c(t), which is also the case for tissue and tumor domains. It is evident that concentration in tumor T1 is higher than in tumor T2 due to the larger diffusion and partitioning coefficients.

FIGURE 6

Figure 6. Concentration evolution in liver. The inlet concentration c(t) at large vessel has a bolus character and generates bolus-type profiles of mean concentrations in capillaries and tissue of the liver and within tumors—reduced with respect to c(t). The lowest concentration is in tumor T2 with the smallest diffusion and partitioning coefficient.

Electrophysiological and Mechanical Model of the Heart Wall

To investigate accuracy of our CSFE model, a small sample of heart wall tissue is selected (Figure 7) following data (Blausen Medical, 2014; Santiago, 2018). From this model, we extract the first layer of muscle cells in myocardium, which is close to the sub-endocardium: the domain that includes the Purkinje fibers.

FIGURE 7

Figure 7. Small domain of heart wall tissue according to Blausen Medical (2014) (Left panel) and first layer of muscle cells close to sub-endocardium with mesh of Purkinje fibers projected on it (Right panel).

According to the image in Figure 7 (right panel), the detailed 2D model is generated (Figure 8A), which consists of the mesh of 1D Purkinje fibers and 25 cells. Dimension of the model is 230 ×150 μm, volume fraction of cells is r_V = 0.71, and area/volume ratio of cell is r_AV = 0.18. Based on the detailed model, we generated the smeared model (Figure 8B) for calculation of the electrical potential and calcium current and concentration (O'Hara et al., 2011). FE nodes of 1D Purkinje fibers are connected with 2D FE nodes of extracellular space domain using connective 1D elements. As shown in Figure 8, it is assumed that the left vertical boundary of the tissue is constrained to displacements in x direction, and lower horizontal boundary is constrained to displacements in y direction. We assumed that muscle fibers have longitudinal direction with respect to the muscle cells, that is, they are aligned to the x direction.

FIGURE 8

Figure 8. (A) The detailed heart wall model with cells and a network of Purkinje fibers. (B) Smeared model with tissue and Purkinje fibers associated to nodes of the CSFEs in a smeared manner.

Data used in the model are: electric conductance (G_i, I = x,y,z) of extracellular space (further called tissue), cells, and in neural fibers is 1,000 [AV⁻¹m⁻¹]; specific membrane conductance (C_m) of Purkinje fibers is 1,000 [S/μm²]; specific membrane capacitance of Purkinje fiber's membrane and cell membrane is 1,000 [AsV⁻¹μm⁻²], and G_m = 0 for the fiber membrane. Diffusion coefficient of Ca²⁺ for tissue and cell is assumed to be 1,000 [μm²/s], while it is assumed that there is no diffusive transport through cell membrane. Material parameters of the muscle mechanical model (Equation 28), used in this example, are: n = 0.4, $C_{50}^n=0.5$, η = 0.2, and σ_max = 100kPa.

The function of the electric potential is taken from Noble (1962): it consists of two identical cycles and is prescribed at inlet nodes of the Purkinje mesh [V_e(t) in Figure 8]. We assumed constant potential inside cells (V_e = −20 mV). Accumulated current density (I_ORd) in cell membrane is calculated according to ORd model (O'Hara et al., 2011) and added to Equation (13) of the FE solution procedure. For these conditions, change of the mean electric potential within tissue is shown in Figure 9A. Results are almost identical for both detailed and smeared models. Currents of ORd model that affect concentration of the Ca²⁺ in myoplasmic compartment are: I_pCa,I_Cab, and I_{NaCa, i}; while currents that affect concentration of Ca²⁺ in the subspace compartment are I_CaL and I_{NaCa, ss}. Mean current density I_Ca for transport of the Ca²⁺ ions can be calculated as (O'Hara et al., 2011):

$\begin{array}{l}{I}_{Ca}=-\left(I_{pCa}+I_{Cab}-2·I_{NaCa,i}\right)\frac{V_{myo}}{V_{myo}+V_{ss}}\\ -\left(I_{CaL}-2·I_{NaCa,ss}\right)\frac{V_{ss}}{V_{myo}+V_{ss}}& (28)\end{array}$

where V_myo and V_ss are volumes of the myoplasmic compartment and subspace compartment, respectively. Change with time of the mean current density obtained by using detailed and smeared model is shown in Figure 9B, where I_ca for the detailed model is calculated as average over all cells.

FIGURE 9

Figure 9. (A) Change of electric potential over time in extracellular space (tissue) domain—detailed and smeared model, with prescribed V_e(t) at inlet nodes of Purkinje fibers and prescribed V_e = −20 V within cells (green). (B) Change of mean current density I_ca [μA/μm²], which affects transport of Ca²⁺ through cell membrane, according to detailed and smeared models.

Mean concentration in cells, ${\left[C{a}^{2+}\right]}_{mean}$, is calculated as average concentration in cells composed of myoplasmic (denoted by index “i”), subspace (“ss”), network SR (“nsr”), and junctional SR (“jsr”) compartments, according to O'Hara et al. (2011)

$\begin{array}{l}{\left[C{a}^{2+}\right]}_{mean}=({\left[C{a}^{2+}\right]}_i·V_{myo}+{\left[C{a}^{2+}\right]}_{ss}·V_{ss}\\ +{\left[C{a}^{2+}\right]}_{nsr}·V_{nsr}+{\left[C{a}^{2+}\right]}_{jsr}·V_{jsr})/V_{cell}& (29)\end{array}$

where V_myo = 0.68V_cell, V_ss = 0.02V_cell, V_nsr = 0.0552V_cell, and V_jsr = 0.048V_cell. Concentrations of Ca²⁺ in each compartment (i, ss, nsr and jsr) of ORd model are calculated according to equations provided in Supplementary of Reference (O'Hara et al., 2011). Change of the mean concentration within cells (cell domain), for both detailed and smeared models, is shown in Figure 10A.

FIGURE 10

Figure 10. (A) Concentration change of Ca²⁺ in cells due to cell membrane currents. (B) Mean contraction (displacement) of the right vertical boundary of heart tissue segment due to Ca²⁺ change in muscle cells.

Muscle contraction occurs from the generation of active stress according to Equation (26), where concentration of calcium is evaluated by our transport models (detailed and smeared). The mechanical response is calculated using the equation of motion (27). Mean contraction (displacement) of the right vertical tissue boundary is shown in Figure 10B. The largest contractions occur at t = 0.9 s and 1.6 s, which are in accordance with the Ca²⁺ concentration within cells.

Effective contraction (modulus of the displacement vector) field of the tissue for the first cycle of action potential function V_e(t) is shown in Figure 11. It can be seen that the largest contraction occurs at t = 0.9 s. Good agreement is observed between the results of the two models.

FIGURE 11

Figure 11. Effective contractions (displacements) according to the detailed model (Left panel) and smeared model (Right panel) for four time points (0.4, 0.9, 1.0, and 1.1 s) of first cycle of action potential function (inlet V_e (t) in Figure 10A).

Electric field potentials within extracellular space for four time points, according to detailed and smeared model, are shown in Figure 12. Potential within cells is kept constant V_e = −20 mV. Agreement between solutions of the two models is noted.

FIGURE 12

Figure 12. Electric potential according to detailed model (Left panel) and extracellular space of smeared model (Right panel) for four time moments (0.4, 0.9, 1.0, and 1.1 s) of the first cycle of action potential function.

Discussion

The main aims of the study were first to demonstrate applicability of our previously published general smeared finite element formulation for the physical fields to large-scale biological problems and second to include electromechanical coupling as the fundamental process for the muscle activation.

Considering the first aim, we note that, by using the today available computational models, software, and hardware, it is practically impossible to calculate drug distribution within an entire organ with different scales, from macro- to nano-, and modeling in detail blood vessels, cells, extracellular space, and even cell interior. On the contrary, our smeared FE methodology offers a platform to achieve solutions for this multiscale problem. We have selected for demonstration a liver of a mouse and generated continuum model with the geometry obtained from images. The model is simple since it consists of 3D continuum finite elements, with the 1D FEs for larger blood vessels as the only addition. The model, for the sake of generality, includes two tumors with different transport properties. The transport parameters data, including hydraulic and diffusive components, are taken according to available references (Gerlowski and Jain, 1986; Sevick and Jain, 1989b, 1991; Less et al., 1991; Mathura et al., 2001; Hall, 2016). The results displayed here show the concentration field, pressure distribution, and fluid velocities in different domains that are quantitatively correct. Accuracy of the solutions can be examined by, first of all, comparison with experiments. Comparison with other available numerical solutions in literature, such as those related to one domain only (e.g., tumor), would be possible with imposing the appropriate boundary conditions in our smeared model. This is not done here, since accuracy of our smeared models was examined in detail in our previous references (Kojic et al., 2017a,b, 2018, 2019; Milosevic et al., 2018b).

Detailed comparison of the smeared models in electrophysiology with the traditional models is given (Kojic et al., 2019). It was found in that analysis that the smeared model is advantageous with respect to the traditional monodomain and bidomain models (e.g., Roth and Wikswo, 1986; Henriquez, 1993; Henriquez et al., 1996; Keener and Panfilov, 1996) (reviewed in Clayton and Panfilov, 2008; Clayton et al., 2011). The defficiency of the traditional models also is that they do not include volumetric fractions of the distinct domains, which obviously must affect the solution. On the other hand, our smeared model relies on the multiple continuum domains coupled by the appropriate connectivity elements. Furthermore, smeared models are simple for application as, for example, in case of the heart: instead to model in detail neural fibers, the entire His-Purkinje system of the heart (Vigmond and Stuyvers, 2016) can be modeled by using the conductivity tensor of the form (19) within continuum 3D finite elements.

With the field of electrical potential within Purkinje network, extracellular and cell interior domain, coupled by the connectivity elements, it is possible to calculate membrane potential that governs the change of the calcium concentration within muscle cells according to Equation (29). It is further straightforward to compute the active stress within muscle cells by employing Equation (26), as one of the commonly used relationships.

Concluding Remarks

The presented smeared finite element formulation offers a platform for multiscale modeling of physical field within a complex medium. Being simple for application, it is suitable for modeling large biological systems, such as liver with tumors used as an example in this study. It is straightforward to implement this approach to other organs, as lung, heart, pancreas, and so on and even brain. The brain has additional complexities related to signal transmission, but the presented model of the Purkinje network and electrical field within tissue shows that a brain model can be developed. The demonstrated electromechanical coupling demonstrates that this methodology can be implemented to muscle modeling in general.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

MK formulated the methodology. MM designed the examples and analyzed results. VS, BM, and VG performed the numerical simulations and prepared results. SN, AZ, NF, and MF provided the data and images for reconstructing the liver and electrophysiology examples.

Funding

This work was supported by the SILICOFCM project that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 777204.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors acknowledge support from the City of Kragujevac, Serbia, Ministry of Education and Science of Serbia (Grants OI 174028 and III 41007), and also from National Cancer Institute under grant U54 CA210187.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2019.00381/full#supplementary-material

References

Baer, S. M., and Rinzel, J. (1991). Propagation of dendritic spikes mediated by excitable spines: a continuum theory. J. Neurophysiol. 65, 874–890. doi: 10.1152/jn.1991.65.4.874

PubMed Abstract | CrossRef Full Text | Google Scholar

Berberoglu, E., Solmaz, H. O., and Göktepe, S. (2014). Computational modeling of coupled cardiac electromechanics incorporating cardiac dysfunctions. Eur. J. Mech. Solids 48, 60–73. doi: 10.1016/j.euromechsol.2014.02.021

CrossRef Full Text | Google Scholar

Blanco, E., Shen, H., and Ferrari, M. (2015). Principles of nanoparticle design for overcoming biological barriers to drug delivery. Nat. Biotechnol. 33, 941–951. doi: 10.1038/nbt.3330

PubMed Abstract | CrossRef Full Text | Google Scholar

Blausen Medical (2014). Medical gallery of Blausen Medical 2014. WikiJ. Med. 1:10. doi: 10.15347/wjm/2014.010

CrossRef Full Text

Canic, S., and Kim, E. H. (2003). Mathematical analysis of the quasilinear e ects in a hyperbolic model blood ow through compliant axi-symmetric vessels. Math. Meth. Appl. Sci. 26, 1161–1186. doi: 10.1002/mma.407

CrossRef Full Text | Google Scholar

Chu, X., Korzekwa, K., Elsby, R., Fenner, K., Galetin, A., Lai, Y., et al. (2013). Intracellular drug concentrations and transporters: measurement, modeling, and implications for the liver. Clin. Pharmacol. Ther. 94, 126–141. doi: 10.1038/clpt.2013.78

PubMed Abstract | CrossRef Full Text | Google Scholar

Cinar, Y., Demir, G., Paç, M., and Cinar, A. B. (1999). Effect of hematocrit on blood pressure via hyperviscosity. Am. J. Hypertens. 12, 739–743. doi: 10.1016/S0895-7061(99)00011-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Clayton, R. H., Bernus, O., Cherry, E. M., Dierckx, H., Fenton, F. H., Mirabella, L., et al. (2011). Models of cardiac tissue electrophysiology: progress,challenges and open questions. Prog. Biophys Mol. Biol. 104, 22–48. doi: 10.1016/j.pbiomolbio.2010.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Clayton, R. H., and Panfilov, A. V. (2008). A guide to modelling cardiac electrical activity in anatomically detailed ventricles. Prog. Biophys Mol. Biol. 96, 19–43. doi: 10.1016/j.pbiomolbio.2007.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Dal, H., Göktepe, S., Kaliske, M., and Kuhl, E. (2013). A fully implicit finite element method for bidomain models of cardiac electromechanics. Comput. Methods Appl. Mech. Eng. 253, 323–336. doi: 10.1016/j.cma.2012.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Apolito, R., Tomaiuolo, G., Taraballi, F., Minardi, S., Kirui, D., Liu, X., et al. (2015). Red blood cells affect the margination of microparticles in synthetic microcapillaries and intravital microcirculation as a function of their size and shape. J. Control Release. 217, 263–272. doi: 10.1016/j.jconrel.2015.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Decker, K. F., Heijman, J., Silva, J. R., Hund, T. J., and Rudy, Y. (2009). Properties and ionic mechanisms of action potential adaptation, restitution, and accommodation in canine epicardium. Am. J. Physiol. Heart Circ. Physiol. 296, H1017–H1026. doi: 10.1152/ajpheart.01216.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrari, M. (2010). Frontiers in cancer nanomedicine: directing mass transport through biological barriers. Trends Biotechnol. 28, 181–188. doi: 10.1016/j.tibtech.2009.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerlowski, L. E., and Jain, R. K. (1986). Microvascular permeability of normal and neoplastic tissues. Microvasc. Res. 31, 288–305. doi: 10.1016/0026-2862(86)90018-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Hall, J. E. (2016). Guyton and Hall Textbook of Medical Physiology, 13th Edn. Philadelphia, PA: Elsevier.

Google Scholar

Henriquez, C. S. (1993). Simulating the electrical behavior of cardiac tissue using the bidomain model. Crit. Rev. Biomed. Engrg. 21, 1–77.

PubMed Abstract | Google Scholar

Henriquez, C. S., Muzikant, A. L., and Smoak, C. K. (1996). Anisotropy, fiber curvature, and bath loading effects on activation in thin and thick cardiac tissue preparations: simulations in a three-dimensional bidomain model. J. Cardiovacs. Eiectrophysiol. 7, 424–444. doi: 10.1111/j.1540-8167.1996.tb00548.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

CrossRef Full Text | Google Scholar

Hunter, P. J., McCulloch, A. D., and ter Keurs, H. E. (1998). Modelling the mechanical properties of cardiac muscle. Prog. Biophys. Mol. Biol. 69, 289–331. doi: 10.1016/S0079-6107(98)00013-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ijiri, T., Ashihara, T., Yamaguchi, T., Takayama, K., Igarashi, T., and Shimada, T. (2008). A procedural method for modeling the Purkinje fibers of the heart. J. Physiol. Sci. 58, 481–486. doi: 10.2170/physiolsci.RP003208

PubMed Abstract | CrossRef Full Text | Google Scholar

Isailovic, V., Kojic, M., Milosevic, M., Filipovic, N., Kojic, N., Ziemys, A., et al. (2014). A computational study of trajectories of micro- and nano-particles with different shapes in flow through small channels. J. Serb. Soc. Comp. Mech. 8, 14–28. doi: 10.5937/jsscm1402014I

CrossRef Full Text | Google Scholar

Jain, R. (1987). Transport of molecules across tumor vasculature. Cancer Metastasis Rev. 6, 559–593. doi: 10.1007/BF00047468

PubMed Abstract | CrossRef Full Text | Google Scholar

Jain, R. (1988). Determinants of tumor blood flow: a review. Cancer Res. 48, 2641–2658.

PubMed Abstract | Google Scholar

Keener, J., and Sneyd, J. (2009). Mathematical Physiology I: Cellular Physiology, II: Systems Physiology, 2nd Edn. New York, NY: Springer.

Google Scholar

Keener, J. P., and Panfilov, A. V. (1996). A biophysical modelfor defibrillation of cardiac tissue. Biophys. J. 71, 1335–1345. doi: 10.1016/S0006-3495(96)79333-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaled, A.-R. A., and Vafai, K. (2003). The role of porous media in modeling flow and heat transfer in biological tissues. Int. J. Heat Mass Transfer 46, 4989–5003. doi: 10.1016/S0017-9310(03)00301-6

CrossRef Full Text | Google Scholar

Kim, N., Cannell, M. B., and Hunter, P. J. (2010). Changes in the calcium current among different transmural regions contributes to action potential heterogeneity in rat heart. Prog. Biophys. Mol. Biol. 103, 28–34. doi: 10.1016/j.pbiomolbio.2010.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Koay, E. J., and Ferrari, M. (2014). Transport oncophysics in silico, in vitro, and in vivo. Preface Phys. Biol. 11:060201. doi: 10.1088/1478-3975/11/6/060201

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojic, M. (2018). Smeared concept as a general methodology in finite element modeling of physical fields and mechanical problems in composite media. J. Serbian Soc. Comput. Mech. 12, 1–16. doi: 10.24874/jsscm.2018.12.02.01

CrossRef Full Text | Google Scholar

Kojic, M., and Bathe, K. J. (2005). Inelastic Analysis of Solids and Structures. Berlin, Heidelberg; New York, NY: Springer.

Google Scholar

Kojic, M., Filipovic, N., Stojanovic, B., and Kojic, N. (2008). Computer Modeling in Bioengineering - Theoretical Background, Examples and Software. Chichester: John Wiley and Sons.

Kojic, M., Mijailovic, S., and Zdravkovic, N. (1998). Modelling of muscle behavior by the finite element method using Hill's three-element model. Int. J. Num. Meth. Engng. 43, 941–953. doi: 10.1002/(SICI)1097-0207(19981115)43:5<941::AID-NME435>3.0.CO;2-3

CrossRef Full Text | Google Scholar

Kojic, M., Milosevic, M., Simic, V., and Ferrari, M. (2014). A 1D pipe finite element with rigid and deformable walls. J. Serb. Soc. Comp. Mech. 8, 38–53. doi: 10.5937/jsscm1402038K

CrossRef Full Text | Google Scholar

Kojic, M., Milosevic, M., Simic, V., Geroski, V., Ziemys, A., and Filipovic, N. (2019). Smeared multiscale finite element model for electrophysiology and ionic transport in biological tissue. Comput Biol. Med. 108, 288–304. doi: 10.1016/j.compbiomed.2019.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojic, M., Milosevic, M., Simic, V., Koay, E. J., Fleming, J. B., Nizzero, S., et al. (2017b). A composite smeared finite element for mass transport in capillary systems and biological tissue. Comp. Meth. Appl. Mech. Engrg. 324, 413–437. doi: 10.1016/j.cma.2017.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojic, M., Milosevic, M., Simic, V., Koay, E. J., Kojic, N., Ziemys, A., et al. (2017a). Extension of the Composite Smeared Finite Element (CSFE) to include lymphatic system in modeling mass transport in capillary systems and biological tissue. J. Serb. Soc. Comp. Mech. 11, 108–120. doi: 10.24874/jsscm.2017.11.02.09

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojic, M., Milosevic, M., Simic, V., Koay, E. J., Kojic, N., Ziemys, A., et al. (2018). Multiscale smeared finite element model for mass transport in biological tissue: from blood vessels to cells and cellular organelles. Comput. Biol. Med. 99, 7–23. doi: 10.1016/j.compbiomed.2018.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojic, M., Slavkovic, R., Zivkovic, M., Grujovic, N., Filipovic, N., and Milosevic, M. (2010). PAK - Finite Element Program for Linear and Nonlinear Analysis. Kragujevac: Univ Kragujevac and R&D Center for Bioengineering.

Kojić, N., Milošević, M., Petrović, D., Isailović, V., Sarioglu, A. F., Haber, D. A., et al. (2015). A computational study of circulating large tumor cells traversing microvessels. Comput. Biol. Med. 63, 187–195. doi: 10.1016/j.compbiomed.2015.05.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Lafortune, P., Arís, R., Vázquez, M., and Houzeaux, G. (2012). Coupled electromechanical model of the heart: parallel finite element formulation. Int. J. Numer. Meth. Biomed. Engng. 28, 72–86. doi: 10.1002/cnm.1494

PubMed Abstract | CrossRef Full Text | Google Scholar

Lauffenburger, D. A., and Horwitz, A. F. (1996). Cell migration: review a physically integrated molecular process. Cell 84, 359–369. doi: 10.1016/S0092-8674(00)81280-5

CrossRef Full Text | Google Scholar

Less, J. R., Skalak, T. C., Sevick, E. M., and Jain, R. K. (1991). Microvascular architecture in a mammary carcinoma: branching patterns and vessel dimensions. Cancer Res. 51, 265–273.

PubMed Abstract | Google Scholar

Ma, P. X., and Zhang, R. (1998). Synthetic nano-scale fibrous extracellular matrix. J. Biomed. Mater. Res. 46, 60–72. doi: 10.1002/(SICI)1097-4636(199907)46:1<60::AID-JBM7>3.0.CO;2-H

PubMed Abstract | CrossRef Full Text | Google Scholar

Mathura, K. R., Vollebregt, K. C., Boer, K., De Graaff, J. C., Ubbink, D. T., and Ince, C. (2001). Comparison of OPS imaging and conventional capillary microscopy to study the human microcirculation. J. Appl. Physiol. 91, 74–78. doi: 10.1152/jappl.2001.91.1.74

PubMed Abstract | CrossRef Full Text | Google Scholar

Mijailovich, S. M., Stojanovic, B., Kojic, M., Liang, A., Wedeen, V. J., and Gilbert, R. J. (2010). Derivation of a finite-element model of lingual deformation during swallowing from the mechanics of mesoscale myofiber tracts obtained by MRI. J. Appl. Physiol. 109, 1500–1514. doi: 10.1152/japplphysiol.00493.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Milosevic, M., Simic, V., Milicevic, B., Koay, E. J., Ferrari, M., Ziemys, A., et al. (2018a). Correction function for accuracy improvement of the composite smeared finite element for diffusive transport in biological tissue systems. Comp. Meth. Appl. Mech. Engrg. 338, 97–116. doi: 10.1016/j.cma.2018.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Milosevic, M., Stojanovic, D., Simic, V., Milicevic, B., Radisavljevic, A., Uskokovic, P., et al. (2018b). A computational model for drug release from PLGA implant. Materials 11:2416. doi: 10.3390/ma11122416

PubMed Abstract | CrossRef Full Text | Google Scholar

Moraru, I. I., Schaff, J. C., Slepchenko, B. M., Blinov, M. L., Morgan, F., Lakshminarayana, A., et al. (2008). Virtual Cell modelling and simulation software environment. IET Syst. Biol. 2, 352–362. doi: 10.1049/iet-syb:20080102

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicholson, C. (2001). Diffusion and related transport mechanisms in brain tissue. Rep. Prog. Phys. 64, 815–884. doi: 10.1088/0034-4885/64/7/202

CrossRef Full Text | Google Scholar

Noble, D. (1962). A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pace-maker potentials. J. Physiolology. 160, 317–352. doi: 10.1113/jphysiol.1962.sp006849

PubMed Abstract | CrossRef Full Text | Google Scholar

Nugent, L. J., and Jain, R. K. (1984). Extravascular diffusion in normal and neoplastic tissues. Cancer Res. 44, 238–244.

PubMed Abstract | Google Scholar

O'Hara, T., Virág, L., Varró, A., and Rudy, Y. (2011). Simulation of the undiseased human cardiac ventricular action potential: Model formulation and experimental validation. PLoS Comput. Biol. 7:e1002061–e1002090. doi: 10.1371/journal.pcbi.1002061

PubMed Abstract | CrossRef Full Text | Google Scholar

Popel, A. S., and Pittman, R. N. (1998). “The microcirculation physiome,” in Handbook of Biomedical Engineering, ed Bronzino, J., (Boca Raton, FL: CRC Press, 911–913.

Google Scholar

Rangamani, P., and Iyengar, R. (2007). Modelling spatio-temporal interactions within the cell. J. Biosci. 32, 157–167. doi: 10.1007/s12038-007-0014-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Reulen, H. J., Graham, R., Spatz, M., and Klatzo, I. (1977). Role of pressure gradients and bulk flow in dynamics of vasogenic brain edema. J. Neurosurg. 46, 24–35. doi: 10.3171/jns.1977.46.1.0024

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, W. G., and Palade, G. E. (1997). Neovasculature induced by vascular endothelial growth factor is fenestrated. Cancer Res. 57, 765–772.

PubMed Abstract | Google Scholar

Rocha, B. M., Kickinger, F., Prassl, A. J., Haase, G., Vigmond, E. J., dos Santos, R.W., et al. (2011). A macro finite element formulation for cardiac electrophysiology simulations using hybrid unstructured grids. IEEE Trans Biomed Eng. 58, 1055–1065. doi: 10.1109/TBME.2010.2064167

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, B. J., and Wikswo, J. P. (1986). A bidomain model for the extracellular potential and magnetic field of cardiac tissue. IEEE Trans. Biomed. Engrg. 33, 467–469. doi: 10.1109/TBME.1986.325804

PubMed Abstract | CrossRef Full Text | Google Scholar

Rushmer, R. F. (1976). Cardiovacular Dynamics, 4th Edn. Philadelphia, PA: W. B. Saunders Co.

Google Scholar

Santiago, A. (2018). Fluid-Electro-Mechanical Model of the Human Heart for Supercomputers. Bacelona: Universitat Politècnica de Catalunya.

PubMed Abstract | Google Scholar

Schaff, J., Fink, C. C., Slepohenko, B., Carson, J. H., and Loew, L. M. (1997). A general computational framework for modeling cellular structure and function. Biophys. J. 73, 1135–1146. doi: 10.1016/S0006-3495(97)78146-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Sevick, E. M., and Jain, R. K. (1989a). Geometric resistance to blood flow in solid tumors perfused ex vivo:Effects of tumor size and perfusion pressure. Cancer Res. 49, 3506–3512.

PubMed Abstract | Google Scholar

Sevick, E. M., and Jain, R. K. (1989b). Viscous resistance to blood flow in solid tumors: effect of hematocrit on intratumor blood viscosity. Cancer Res. 49, 3513–3519.

PubMed Abstract | Google Scholar

Sevick, E. M., and Jain, R. K. (1991). Effect of red blood cell rigidity on tumor blood flow: increase in viscous resistance during hyperglycemia. Cancer Res. 51, 2727–2730.

PubMed Abstract | Google Scholar

Shi, Y., Lawford, P., and Hose, R. (2011). Review of zero-D and 1-D models of blood flow in the cardiovascular system. Biomed. Eng. 10:33. doi: 10.1186/1475-925X-10-33

PubMed Abstract | CrossRef Full Text

Skinner, S. A., Tutton, P. J., and O'Brien, P. E. (1990). Microvascular architecture of experimental colon tumors in the rat. Cancer Res. 50, 2411–2417.

PubMed Abstract | Google Scholar

Slepchenko, B. M., Schaff, J. C., Macara, I., and Loew, L. M. (2003). Quantitative cell biology with the Virtual Cell. Trends Cell Biol. 13, 570–576. doi: 10.1016/j.tcb.2003.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, N. P., Pullan, A. J., and Hunter, P. J. (2002). An anatomically based model of transient coronary blood flow in the heart. SIAM J. Appl. Math. 62, 990–1018. doi: 10.1137/S0036139999355199

CrossRef Full Text | Google Scholar

Stojanovic, B., Kojic, M., Rosic, M., Tsui, C. P., and Tang, C. Y. (2007). An extension of Hill's three-component model to include different fiber types in finite element modelling of muscle. Int. J. Num. Meth. Eng. 71, 801–817. doi: 10.1002/nme.1963

CrossRef Full Text | Google Scholar

Sun, H., and Pang, K. S. (2008). Transport, and metabolism of solutes in Caco-2 cell monolayers: a theoretical study. Drug Metabol. Dispos. 36, 102–123. doi: 10.1124/dmd.107.015321

PubMed Abstract | CrossRef Full Text | Google Scholar

Swabb, E. A., Wei, J., and Gullino, P. M. (1974). Diffusion and convection in normal and neoplastic tissues. Cancer Res. 34, 2814–2822.

PubMed Abstract | Google Scholar

Trapp, S., and Horobin, R. W. (2005). A predictive model for the selective accumulation of chemicals in tumor cells. Eur. Biophys. J. 34, 959–966. doi: 10.1007/s00249-005-0472-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigmond, E. J., and Stuyvers, B. D. (2016). Modeling our understanding of the His-Purkinje system. Prog. Biophys. Mol. Biol. 120, 179–188. doi: 10.1016/j.pbiomolbio.2015.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Winslow, R. L. (1992). Theoretical Foundations of Neural Modeling, BME 580.681. Baltimore: The John Hopkins Univ. Sch. Medicine and Whiting Sch Engrg.

Zagorchev, L., Oses, P., Zhuang, Z. W., Moodie, K., Mulligan-Kehoe, M. J, Simons, M., et al. (2010). Micro computed tomography for vascular exploration. J Angiogenes Res. 2:7. doi: 10.1186/2040-2384-2-7

PubMed Abstract | CrossRef Full Text | Google Scholar