Tissue perfusion model

[menon-karthik]

I am moving this issue over verbatim from the Fortran svFSI repository

I am developing a model for the perfusion of blood within a tissue volume. The model is a single-compartment Darcy model, which can be converted into a Helmholtz equation.

Details of the model in this paper: Papamanolis, L., et al. (2021). Myocardial Perfusion Simulation for Coronary Artery Disease: A Coupled Patient-Specific Multiscale Model. Annals of Biomedical Engineering, 49(5), 1432–1447. https://doi.org/10.1007/s10439-020-02681-z

Details of implementation:

The solver for the Helmholtz equation is being developed as a modification to the HEATS diffusion solver.
The tissue volume is divided into perfusion regions, each with a perfusion pressure. Other important parameters are permeability and conductance parameters.
The nodes corresponding to each perfusion region, and the corresponding perfusion pressure, depend on the vascular network that is perfusing the tissue volume. So this is user-specific, but the interface with svFSI is a list of nodes that specify each territory.
This is generally solved as a steady problem, which can easily be accomplished by setting density = 0.
The idea is to couple this with a flow simulation so that the flow pressures drive the perfusion in the tissue volume. The aim is to couple it with 3D and 0D flow solvers. Currently, the interface between the perfusion model and the flow solver will be a "hand-off" using file I/O (subject to improvement if required).

[zasexton]

I'll be looking into this issue this week to see if we can finish moving the perfusion model over and apply it to a few test cases. @ktbolt would you mind assigning this issue to me and @menon-karthik ?

[ktbolt]

@zasexton Can't you assign it to yourself?

[zasexton]

@ktbolt I don't think so...

[zasexton]

To better help me understand how to implement the perfusion equation in svFSI in a consistent and true manner I'll be including my brief formulation notes here. Once we agree on the formulation, I'll type up the full derivation for future students/developers becoming familiar with new equations in svFSI.

Ultimately I would like to implement the perfusion equations in an unsteady manner but so far I have the following...

Strong Form of the Unsteady Single-Compartment Darcy Model

$$ \dfrac{\partial \boldsymbol{u} }{\partial t} + \boldsymbol{u} + K \nabla p = 0$$

$$\nabla \cdot \boldsymbol{u} = \beta_{\text{source}}(p_{\text{source}} - p) - \beta_{\text{sink}}(p-p_{\text{sink}})$$

We further assume that permeability is homogeneous and isotropic (i.e. $K = k$ thus it can be treated as a constant value). $\boldsymbol{u}$ is the Darcy volume flux vector. $p$ is the pressure. $p_{\text{source}}$ and $p_{\text{sink}}$ are the pressure sources and sinks applied as Dirichlet boundary conditions, respectively. $\beta_{\text{source}}$ and $\beta_{\text{sink}}$ are the admittance parameters relating the pressure gradient to the volume flux.

Weak Form of the Unsteady Single-Compartment Darcy Model

Let $\mathcal{S}$ and $\mathcal{Q}$ denote the trial and variational spaces consisting of real-valued functions of sufficient smoothness and meet the requirements imposed from the strong form. Given $\beta_{\text{source}}, p_{\text{source}}, \beta_{\text{sink}}, p_{\text{sink}}, k \rightarrow \mathbb{R}$, find $p \in \mathcal{S} $ and $q \in \mathcal{Q}$ such that;

$$\int_{\Omega} \nabla \cdot \dfrac{\partial \boldsymbol{u}}{\partial t} qd\Omega -\int_{\Omega}k \nabla p \nabla q d\Omega + (\beta_{\text{source}} + \beta_{\text{sink}}) \int_{\Omega} pq d\Omega = - \int_{\Omega} (\beta_{\text{source}} p_{\text{source}} + \beta_{\text{sink}} p_{\text{sink}}) d\Omega - \int_{\Gamma} k (\nabla p \cdot \boldsymbol{n})q d\Gamma $$

Code implementation for the Steady 2D Darcy Problem

If we neglect the time components of the previous weak form and consider only the 2D case, we may implement the full-discrete, finite element formulation of the previous weak problem. I denote $\beta_{\text{source}}$ as b0 and $\beta_{\text{sink}}$ as b1. The $p_{\text{source}}$ and $p_{\text{sink}}$ components are represented by local_source and local_sink, respectively, for the local residual lR and local stiffness lK contributions to the global residual and stiffness matrices, respectively. We neglect the local contribution to the acceleration matrix al due to our assumption of steady-state.

        double Td = -source;
        Vector<double> Tx(nsd);


        for (int a = 0; a < eNoN; a++) {
            Td = Td - local_source(a)*N(a)*b0 - local_sink(a)*N(a)*b1;
            Td = Td - (b0 + b1)*yl(i,a);
            Tx(0) = Tx(0) + Nx(0, a) * yl(i, a);
            Tx(1) = Tx(1) + Nx(1, a) * yl(i, a);
        }
        
        for (int a = 0; a < eNoN; a++) {
            lR(0, a) = lR(0, a) + w * (N(a) * Td  - k * (Nx(0, a)*Tx(0)
                    + Nx(1, a) * Tx(1)));
            for (int b = 0; b < eNoN; b++) {
                lK(0, a, b) = lK(0, a, b) + wl * (N(a) * N(b) * (b0 + b1) -
                        k * (Nx(0, a) * Nx(0, b) + Nx(1, a) * Nx(1, b)));
            }
        }

[zasexton]

So I've been trying to add a new perfusion equation to svFSIplus. For now it is essentially a duplicate of the heats.cpp and heats.h files which have names darcy.cpp and darcy.h, respectively. My goal is to initially reproduce the result solutions of the line source diffusion test case using the EquationType::phys_darcy so that I can be sure that the data handling is correct before changing any of the physics.

When I have done this so far I have only gotten zero solutions:

Darcy Solution:

Expected Solution (from Heats.cpp):

Obviously there is an issue with the LR or LK values that fails to update from zero but I am unsure how to resolve this. So far I have updated CMakeLists.txt, ComMod.h, consts.cpp, consts.h, darcy.cpp, darcy.h, eq_assem.cpp, initialize.cpp, Parameters.cpp, Parameters.h, read_files.cpp, set_equation_dof.h, and set_equation_props.h to reflect only the minimal changes needed to build svFSIplus with a new equation that essentially is a duplicate of the heats implementation. @ktbolt Are there potentially other files that I am forgetting to update in order to carry these changes through?

[zasexton]

Found the error in the darcy.cpp file with an incorrect bracket closure. This led to incorrect assembly of lR and lK assembly among the elements. I have commited a new change to resolve this and the expected solution is obtained from the darcy equation and matches the heatS equation.

Darcy Solution:

Expected Solution: (HeatS Equation)

[menon-karthik]

To better help me understand how to implement the perfusion equation in svFSI in a consistent and true manner I'll be including my brief formulation notes here. Once we agree on the formulation, I'll type up the full derivation for future students/developers becoming familiar with new equations in svFSI.

Ultimately I would like to implement the perfusion equations in an unsteady manner but so far I have the following...

Strong Form of the Unsteady Single-Compartment Darcy Model

∂u∂t+u+K∇p=0

∇⋅u=βsource(psource−p)−βsink(p−psink)

We further assume that permeability is homogeneous and isotropic (i.e. K=k thus it can be treated as a constant value). u is the Darcy volume flux vector. p is the pressure. psource and psink are the pressure sources and sinks applied as Dirichlet boundary conditions, respectively. βsource and βsink are the admittance parameters relating the pressure gradient to the volume flux.

Weak Form of the Unsteady Single-Compartment Darcy Model

Let S and Q denote the trial and variational spaces consisting of real-valued functions of sufficient smoothness and meet the requirements imposed from the strong form. Given βsource,psource,βsink,psink,k→R, find p∈S and q∈Q such that;

∫Ω∇⋅∂u∂tqdΩ−∫Ωk∇p∇qdΩ+(βsource+βsink)∫ΩpqdΩ=−∫Ω(βsourcepsource+βsinkpsink)dΩ−∫Γk(∇p⋅n)qdΓ

Code implementation for the Steady 2D Darcy Problem

If we neglect the time components of the previous weak form and consider only the 2D case, we may implement the full-discrete, finite element formulation of the previous weak problem. I denote βsource as b0 and βsink as b1. The psource and psink components are represented by local_source and local_sink, respectively, for the local residual lR and local stiffness lK contributions to the global residual and stiffness matrices, respectively. We neglect the local contribution to the acceleration matrix al due to our assumption of steady-state.

        double Td = -source;
        Vector<double> Tx(nsd);


        for (int a = 0; a < eNoN; a++) {
            Td = Td - local_source(a)*N(a)*b0 - local_sink(a)*N(a)*b1;
            Td = Td - (b0 + b1)*yl(i,a);
            Tx(0) = Tx(0) + Nx(0, a) * yl(i, a);
            Tx(1) = Tx(1) + Nx(1, a) * yl(i, a);
        }
        
        for (int a = 0; a < eNoN; a++) {
            lR(0, a) = lR(0, a) + w * (N(a) * Td  - k * (Nx(0, a)*Tx(0)
                    + Nx(1, a) * Tx(1)));
            for (int b = 0; b < eNoN; b++) {
                lK(0, a, b) = lK(0, a, b) + wl * (N(a) * N(b) * (b0 + b1) -
                        k * (Nx(0, a) * Nx(0, b) + Nx(1, a) * Nx(1, b)));
            }
        }

Just to expand on this a little, I had implemented a version of this in the svFSI Fortran code. I re-implemented it in svFSIplus a couple of months ago, but could not get svFSIplus to reproduce the results from the Fortran code. I did not investigate very closely at the time but that is ongoing.

About the model itself, the current plan is to implement this for a steady Darcy problem. This is what is relevant to most applications. We will explore the unsteady formulation at a later time. So what we are going to be solving is:

$$\vec{u} + K\nabla P = 0$$

$$ \nabla \cdot \vec{u} = s(P)$$

We can use the second equation to eliminate $\vec{u}$, and so what we will really solve is

$$K \nabla^2 P = -s(P)$$

This looks like a steady diffusion equation. But since $s(P)$ depends on $P$, in reality it turns out to be closer to something of the form:

$$K \nabla^2 P + \alpha P= \gamma$$

Since it is similar to a diffusion equation, we are going to implement this by starting with the heats model and adding in the extra terms.

[vvedula22]

Currently, none of our equations are written or solved for steady state. In general, we solve for the steady state problem by assuming a pseudo time derivative and solve the resulting system as an unsteady problem that converges to a steady state (i.e., variables don't change with time up to a suitable tolerance). Since you are planning to have an unsteady formulation in the future, it is best to implement an unsteady model now and consider the steady state problem as a special case.

If you agree, IMO, this problem is best solved using a modified Navier-Stokes equation (most general form) or the Stokes equation (neglecting inertial effects), since you have a modified form of the continuity equation (with a pressure-dependent source) term coupled to the modified momentum equation (with a velocity-dependent source term). You might also have to deal with stabilization for pressure-velocity coupling. Therefore, I think you can easily implement this by modifying the FLUID or STOKES codes by adding suitable domain-dependent parameters. These parameters will trigger the corresponding source terms in the momentum and continuity equations to convert them into the Navier-Stokes-Brinkman equations (I believe this model is the general model for Darcy-type perfusion).

[zasexton]

@vvedula22 this would definitely be good to implement and will likely be a better solution to the family of flow problems likely to be encountered from tissue engineering applications. Since I do not currently have the level of familiarity with svFSI/svFSIplus architectures I will likely plan to continue implementing the simpler perfusion model for now.

This may also be advantageous since this model has recently been appreciated for heart perfusion modeling (@menon-karthik
's original post); thus if the proposed Navier-Stokes-Brinkman (NSB) equations are more accurate to the physics we will have the simple alternative model to test against (even though as you mention steady Darcy flow is surely a special case from the NSB equations).