Isogeometric Analysis for Fluid Shear Stress in Cancer Cells

The microenvironment of the tumor is a key factor regulating tumor cell invasion and metastasis. The effects of physical factors in tumorigenesis is unclear. Shear stress, induced by liquid flow, plays a key role in proliferation, apoptosis, invasion, and metastasis of tumor cells. The mathematical models have the potential to elucidate the metastatic behavior of the cells’ membrane exposed to these microenvironment forces. Due to the shape configuration of the cancer cells, Non-uniform Rational B-splines (NURBS) lines are very adequate to define its geometric model. The Isogeometric Analysis allows a simplified transition of exact CAD models into the analysis avoiding the geometrical discontinuities of the traditional Galerkin traditional techniques. In this work, we use an isogeometric analysis to model the fluid-generated forces that tumor cells are exposed to in the vascular and tumor microenvironments, in the metastatic process. Using information provided by experimental tests in vitro, we present a suite of numerical experiments which indicate, for standard configurations, the metastatic behavior of cells exposed to such forces. The focus of this paper is strictly on geometrical sensitivities to the shear stress’ exhibition for the cell membrane, this being its innovation.


Introduction
The formation of a secondary tumor at a site distant from the the primary tumor is known as metastasis by a cancerous tumor. To initiate the metastatic spread of cancer through the bloodstream, tumor cells must transit through microenvironments of dramatically varying physical forces. Cancer cells are able to migrate through the stroma, intravasate through the endothelium into the blood or lymphatic vessels, to flow within the vessels, to extravasate from the vessel through the endothelium and colonize in tissue at a secondary site [1].
In soft tissues, cancer cells are exposed to mechanical forces due to fluid shear stress, hydrostatic pressure, tension and compression forces. Fluid shear stress is one of the most important forces that cells are exposed to, and its effects on blood cells, endothelial cells, smooth muscle cells, and others have been extensively studied. However, much less is known about fluid shear stress effects on tumor cells. Cancer cells experience two main kinds of fluid shear stress: stresses generated by blood flow in the vascular microenvironment, and those generated by interstitial flows in the tumor microenvironment [2]. Stresses generated by interstitial and blood flows could contribute to the metastatic process by enhancing tumor cell invasion and circulating tumor cell adhesion to blood vessels, respectively. However, it is difficult to predict tumor cell behavior to such forces and it is difficult to experimentally measure such flows in the tumor microenvironment [3].
In this work, we apply Mathematical and Mechanical processes to analyze the metastatic behavior of the cells' membrane exposed to microenvironement forces. Providing an Isogeometric Analysis (IGA) [4] computational method to model and predict how cancer cells respond to such forces, we allow for new insight and new decision tools for medical problems.
The biological fluid study at microenvironments, in vivo or in vitro, is a recent topic following the microfluidic technology and mechanical methodology methods developments [5][6][7][8], with the ethical and budget restrictions, increasing the evidence that fluid shear stress is an essential factor affecting fluid mechanics [1].
The current work serves as an introduction to this line of research on mathematical modeling to help us understand the omics data produced by experimental techniques and to bridge the gap between the developments of technologies and systemic modelling of the biological process in cancer research.
Although it is accepted that the effectiveness of IGA methods is well-established for problems with complex geometries, its effect in biologic models has not been extensively studied. Indeed, at present, general research in this subject is still in this infancy for the case studies presented [9,10].
First, we introduce and describe the methodology. We then provide a mathematical and numerical analysis overview for the interstitial flow governing model and in order to estimate the fluid shear stress on cells in tissues.
Standard shape for cancer cells is very irregular, for this reason Non-uniform Rational B-splines (NURBS) lines are very adequate to define its geometric model. The Isogeometric Analysis allows for simplified transition of exact CAD models into the analysis, avoiding the geometrical discontinuities of the traditional Galerkin traditional techniques. For the exhibit models, NURBS provides a high-quality geometric model, quite analogous to a physical model, and also defines the basis of the discrete space in which the partial differential equation solution is approximated with great accuracy per degree of freedom, as referred to in Section 3.
We conclude with simulations to predict the effect of fluid shear stress on cancer cells for several scenarios of microenvironment tumor implantation.
In this paper, we focus on the effects of shear stress, induced by liquid flow, on apoptosis, invasion, and metastasis of tumor cells, following the theoretical [11] and numerical [12] work for a Stokes flow.

Model Geometry
Tumor cells have a large variety in shape and size. Therefore, several different geometries will be considered in this study, always considering these as standard cells. We will consider bidimensional projections (as shown in Figure 24) of the model shown in Figure 1 to inspire the geometrical model for the presented cases.  [7]. Bidimensional projections of this model will be considered for the study cases.
For these geometrical models, Non-uniform Rational B-splines (NURBS) lines are very adequate to define it due to simplified transition of exact CAD models into the analysis, avoiding the geometrical discontinuities of the traditional Galerkin traditional techniques.

Governing Equations
Average equations describing viscous flow through a porous region are of great theoretical and practical interest. At the microscale level, the Stokes equations apply and provide a complete description of the entire flow field. However, Darcy's law, formally derived by performing appropriate volume averages of Stokes equations, is applicable. The qualitative difference between these two flow descriptions motivate Brinkman to suggest a general equation that interpolates between the Stokes equation and Darcy's law [13]. We introduce and describe the mathematical model in the next section.

The Darcy-Brinkman Equation
Darcy's law is a phenomenologically derived constitutive equation that describes the flow of a fluid through a porous medium [9,10]. This law, as an expression of conservation of momentum, was determined experimentally by Darcy. It has since been derived from the Navier-Stokes equations via homogenization. It is analogous to Fourier's law in the field of heat conduction.
where k is the permeability of the medium, ∇p is the pressure gradient vector, ν is the viscosity of the fluid and u is the fluid's average velocity through a plane region represented by Ω. Region Ω is supposed regular enough to ensure the later theoretical regularity for u and p.
To account for interstitial flows between boundaries, Brinkman has developed a second-order term, taking into account no-slip boundary conditions cells' membrane ( Figure 2). For governing interstitial flow between boundaries we will consider the Darcy-Brinkman equation (2) The mass balance equation for a steady state incompressible fluid is that the divergence of the fluid is zero To apply the boundary conditions we decompose the boundary of the region Ω in two non-overlaping regions, as shown in Figure 3, ∂Ω = Γ T ∪ Γ I , the boundary of the tumor region and the boundary of the interstitial region considered, respectively. We consider homogeneous Dirichlet boundary condition for the velocity over the cells membrane and non-homogeneous on the interstitial region boundaries, i.e., u = (u 0 x , u 0 y ) and parameters ν and k given by Table 1, for the cases under study. To facilitate understanding, we write the component-wise formulation in two dimensions. Let u = (u x , u y ). Then Equations (2) and (3) consists of three equations We are interested to find u = (u x , u y ) and p solution of the problem defined by Equation (4) and the Dirichlet boundary conditions showed at Table 1.
With the solution, u = (u x , u y ) we evaluate the fluid shear stress. For bidimensional Newtonian fluid the Stress is written [14] and the Shear Stress is described by

Isogeometric Discretization
In order to represent complex shapes, the use of polynomials or rational segments may often be inadequate or imprecise. On the other hand, B-spline and NURBS functions enjoy some major advantages that make them extremely convenient for complex geometrical representations.
The main idea behind the isogeometric approach [4] is to discretize the unknowns of the problem with the same set of basis functions that CAD employs for the construction of geometries. Let p be the prescribed degree and n control points, we define by the knots vector, with t i ∈ [0, 1], i ∈ {1, · · · , n + p + 1}. Cox-De Boor's formula [15] defines n one-dimensional B-spline basis functions recursively as As referred to in [12] the support of a B-spline of degree p is always p + 1 knot spans and, as a consequence, each p-th degree function has p − 1 continuous derivatives across the element boundaries, or across the knots, if they are not repeated. Repetition of knots can be exploited to prescribe the regularity.
NURBS of degree p are defined as rational B-splines, associating to each B-spline function a weight w i with w i called the weight parameter. Geometries in the projective space may be described by using the concept of homogeneous coordinates, which are frequently denoted as weights w i . A weighted polynomial B-spline geometry of IR d+1 is obtained by first multiplying its control point data with the homogeneous coordinates. For values w i > 1, the object moves toward the control polygon, whereas for weights smaller than one, the influence of the control point on the geometry decreases. Control points with w i = 0 do not affect the geometric object at all. If all w i are equal to one, the NURBS basis simplifies to the polynomial B-spline basis [15] and allows for the partition of unity property. As in [12] we define bidimensional B-splines and NURBS using a tensor product approach. Considering Ξ = Ξ 1 × Ξ 2 the knot vectors, p the degrees and n the number of basis functions the bivariate B-spline is given by where t = (t 1 , t 2 ), p = (p 1 , p 2 ), n = (n 1 , n 2 ) and i = (i 1 , i 2 ) is a multi-index in the set i 1 ∈ {1, · · · , n 1 } and i 2 ∈ {1, · · · , n 2 }. It is straightforward to notice that there is a parametric Cartesian mesh T h associated with Ξ. The knot vectors partitioning the parametric domain [0, 1] 2 into parallelograms. For each element Q ∈ T h we associate a parametric mesh size h Q = diam(Q), and h = max h Q , Q ∈ T h .
In the following, we refer to the basis functions indicating the global index, and we will denote by S h (Ξ) the bivariate B-spline space spanned by the basis functions N i , 1 ≤ i ≤ n. For convenience we also use the notation S p 1 ,p 2 α 1 ,α 2 (Ξ) to designate the associated space of splines of order p 1 in the x direction, p 2 in the y direction, and smoothness α 1 and α 2 respectively. Notice that α i is determined by the knot vector Ξ i , spanned by the basis functions of degree p d and regularity α d for each direction d ∈ {1, 2}.
A NURBS curve is defined by a set of control points P which act as weights for the linear combination of the basis functions, giving the mapping to the physical space. In particular, given n one-dimensional basis functions N i,p and n control points P i ∈ IR 2 , i ∈ {1, · · · , n} a curve parameterization is given by: The control points define the so-called control mesh but this does not, in general, conform to the actual geometry (cf Figure 5b). In many real-world applications, the computational domain may be too complicated to be represented by a single NURBS mapping from the reference domain to the physical space. This could be due to topological reasons (or to the presence of different materials). In these cases it is common practice to resort to the so-called multipatch approach [4,12] (cf Figure 4). Here, the physical domain is split into simpler subdomains Ω k such as Ω = ∪ k Ω k and Ω i ∩ Ω j = ∅, or a point, or an edge, for i = j.

Isogeometric Conforming Spaces
Analogously to classical Finite Element Methods, IGA is based on a Galerkin approach: the equations are written in their variational formulation, and the solution is sought in a finite-dimensional space with the correct approximation properties. In IGA the basis function space is inherited from the one used to parametrize the geometry.
Let us now consider a domain Ω ⊂ V that can be exactly parametrized with a mapping F For a multipatch approach we will consider a map for each Ω k . Then, the discrete space in the physical domain is defined applying the isoparametric concept as where each Π has the same knots Ξ but the multiplicity has been increased by one, that means the velocity components' space have the same continuity as the pressure space [16]. We recall the important property of B-spline: at a knot of multiplicity m, basis function N i,p is C p−m = C α continuous. Proposition 1. With the notation and assumptions above, the space is subspace of Q h .

Proof of Proposition 1. For
The subspace V 0 h ⊂ V h h is the space of discrete functions that vanish on the boundary of Ω.

The Variational Formulation Discretization
We consider a classic mixed variational discretization of problem (4) in primitive variables [16,17], in which an approximation ( u h , p h ) to the exact solution ( u, p) of (4) is obtained by solving the problem where (V h , M h ) are couples of finite-dimensional spaces parameterized with h, as introduced above, and with the bilinear forms a and b defined as The conditions for the well posedness of a saddle point system is known as inf-sup conditions or Ladyzhenskaya-Babuška-Breezi (LBB) condition: with C independent of the discritization parameter h, · L 2 · H 1 are the classic norms defined over the spaces L 2 (Ω) and H 1 (Ω), respectively, as in [11]. The subspace choice M h ensures the discrete velocity u h solution of (16) is divergence-free.
In the context of IGA, we will use spline-based spaces V 0 h , M h , which satisfy the inf-sup condition (17). The elements considered can be seen as spline generalization of well-known finite elements spaces, namely Taylor-Hood elements [11]. Thanks to the high interelement regularity, which is the main feature of splines, the proposed discretizations are conforming, i.e., produce globally continuous, discrete velocities. Moreover, the spline generalization of Taylor-Hood elements also enjoys property [18] and thus provides divergence-free discrete solutions (15).

Numerical Results
In this section, we present some numerical experiments to predict the effect of fluid shear stress on cancer cells localized in the interstitial region to assess its metastatic behavior. For all cases, we use experimental data obtained in the works of [1,[5][6][7][8].
We use GeoPDEs, an open source and free package for isogeometric analysis in Matlab [12,19], to achieved the numerical implementation of the discretized problem (16).
As a peculiarity of IGA is to allow for high degree and high regularity discretization spaces, most of the tests below are performed for degree p = 5 and regularity α = 4.
We use multipatch geometries ( Figure 4) and these discretizations are C 0 between the patches. When applied to incompressible flows, these discretizations produce pointwise divergence-free velocity fields and hence exactly satisfy mass conservation. We enforce the Dirichlet boundary conditions weakly by Nitsche's method, allowing this method to default to a compatible discretization of Darcy flow in the limit of vanishing viscosity [12].
Using the data from Table 1 we perform different experiments with different choices of region configuration and tumor size.

Case 1
Let us begin by considering a simple circular case which will be used as a pattern for the others ones. We consider Ω = [−10, 10] × [−10, 10] and a centred tumor region with a circular boundary, with radius r = 2(µ).
For this case we present the physical initial mesh, with control points, at Figure 5, an h-refinement (by multiple knot insertion) sequence, Figures 6 and 7, and the correspondent results for velocity and shear stress.  With this discretization the velocity magnitude and stream lines obtained is represented in Figure 8. The flow-generate shear stress is between −7 × 10 −2 and 7 × 10 −2 dyn/cm 2 (we recall 1 dyn/cm 2 = 0.1 Pa) over the whole domain, as shown in Figure 9, and maximum is achieved on the cell surface, this means on the membrane as expected. Next, for this case, we will repeat the simulation with a h-refinement on the domain. First we consider at patches interfaces the new knots array Ξ h 1 = [0.8, 0.95] .
With this discretization the flow-generate shear stress is between −1.5 and 1.5 dyn/cm 2 over the whole domain, as shown in Figure 10, and velocity magnitude is between 0 and 2 µs −1 , Figure 11.  With this discretization the flow-generate shear stress is between −1.5 and 1.5 dyn/cm 2 over the whole domain, as shown in Figure 12, and velocity magnitude is between 0 and 2 µs −1 , Figure 13, as in the previous refinement case.  These numerical results show that mesh refinement does not affect results for T h 2 . At this point, the model and its results are independent of the mesh.
Finally, we point to the velocity magnitude computed at any point P ∈ Ω using discrete 2 norm

Case 2
Case 2 is similarly to case 1, we consider Ω = [−10, 10] × [−10, 10] and a centred tumor region with a circular irregular boundary, as shown at Figure 14. Whereas the irregular circular line is obtained from perturbations on the line of the pattern case 1, we quantify the boundary irregularities by comparison between the cells' perimeter, introducing the parameter I (the irregularity parameter indexed with the case study number) defined as with d 1 = 4π = 12.5664 and d 2 = 13.3937. For Case 2 we obtain I 2 = 0.0658. With this cell configuration, I 2 = 0.0658, the flow-generate shear stress is between −1.7 and 1.7 dyn/cm 2 over the whole domain, as shown in Figure 15, and velocity magnitude is between 0 and 2.2 µs −1 , Figure 16.

12.5664
= 0.1014, as shown in Figure 17. With this cell configuration, I 3 = 0.1014, the flow-generate shear stress is between −2.2 and 2.2 dyn/cm 2 over the whole domain, as shown in Figure 18, and velocity magnitude is between 0 and 2.5 µs −1 , Figure 19.   With this cell configuration, I 4 = 0.613, the flow-generate shear stress is between −3 and 3 dyn/cm 2 over the whole domain, as shown in Figure 21, and velocity magnitude is between 0 and 2.9 µs −1 , Figure 22.  With the results from 1 to 4 we can define an evolution graphic in Figure 23 for shear stress maximum, the blue line, and velocity magnitude maximum, the red line. Regarding cell membrane resistance, at microenvironment context this evolution is quite important. We observed a linear dependence between the irregularity parameter and shear stress. The Figure 23 shows the dependence of shear stress on velocity and the irregularity parameter, given that the two lines have different growth rates. To conclude we present two more cases quite near to reality. Inspired by Figure 1 we define the outline shown in Figure 24. The following cases are determined from the outline 3D by cut planes for case 5, the green one, and case 6, the blue one.

Case 5
With this cell configuration, domain in Figure 25 (corresponding to the green section in Figure 24) we get I 5 = 0.613, the flow-generate shear stress is between −2.8 and 2.5 dyn/cm 2 over the whole domain, as shown in Figure 26, and velocity magnitude is between 0 and 2.8 µs −1 , Figure 27. For this value of I 5 we obtain the expected from comparison with the above cases. We point to the geometric asymmetric effect on the maximum value of shear stress.

Case 6
Finally we consider for the cell configuration the blue section of the model in Figure 24, as shown in Figure 28. The flow-generate shear stress is between −4.1 and 3.3 dyn/cm 2 over the whole domain, as shown in Figure 29, and velocity magnitude is between 0 and 1.9 µs −1 , Figure 30.
With this result, we can see that the zone with the highest shear stress value is over the "body" of the cell and not in its branches. In fact, it is the branches we should focus on because they promote the displacement or duplication of the cell.

Conclusions
From this preliminary study, it is possible to establish the adequacy of the isogeometric analysis and NURBS to model irregular tumor cells and evaluate the shear forces at tumor microenvironments.
With some numerical examples, using different cell configurations, we have identified conditions that allow the increase or decrease of the fluid shear stress, which could contribute to the metastatic process by enhancing tumor cell invasion and circulating tumor cells. The cancer cell invasive potential is significantly reduced, as much as 92%, upon exposure to 0.55 dyn/cm 2 fluid shear stress [21]. This work contributed to establish a relationship between a quantification of cell membrane irregularity and the maximum value of shear stress evaluated on the cell membrane, by the effect of interstitial fluid in the microenvironment. In follow-up works, we must understand how the cell reacts to these forces.
Mathematical models of fluid shear stress effects coupled with in vitro and in vivo experimental validation, may better predict cell behavior in such dynamic microenvironments, and potentially provide novel approaches for the prevention of metastasis.
A particular feature of cancer modeling revolves around the idea that a tumors' size and shape changes over time and its resistance to the fluid shear stress is also affected by the presence of other substances in the interstitial region meaning, for this reason, is the motivation for further work.