User Guide

This guide provides a brief overview of the package and its usage. For more detailed information, refer to the API documentation.

Installation

The recommended way to install HemTracer is to clone the GitHub repository and install the package using a virtual environment.

git clone https://github.com/nicodirkes/HemTracer.git
cd HemTracer
python3.12 -m venv .env
source .env/bin/activate
python3.12 -m pip install .

The virtual environment has to be activated whenever you want to use HemTracer. Note that the package requires Python 3.12 due to typing functionality. In addition, the package requires the following dependencies:

• numpy

• scipy

• vtk

These dependencies are automatically installed when installing the package using pip. The versions used during development and testing are listed in requirements.txt in the repository.

Usage

The package is meant to be used as part of a Python script. After installation, it can be imported as follows:

import hemtracer as ht

Importing CFD data

The first step is to read in the CFD data. This is handled by the hemtracer.EulerianFlowField class. The constructor takes the path to the VTK file as an argument:

flow_field = ht.EulerianFlowField('path/to/vtk/file.vtk')

If the CFD data represents an MRF simulation (as should be the case for most centrifugal blood pumps), the velocities have to be given all in the absolute frame. The package will then perform the necessary transformations using the hemtracer.EulerianFlowField.transform_mrf_data() method:

flow_field.transform_mrf_data(rf_name='ref_frame', rf_rot=2, omega=[0, 100, 0], x0=[0,0,0])

Refer to the documentation of the hemtracer.EulerianFlowField.transform_mrf_data() method for more details on the arguments. Note that this step is required for all MRF simulations.

Computing pathlines

The next step is to compute pathlines. This is handled by the hemtracer.PathlineTracker class. The constructor takes the flow field as an argument:

pathline_tracker = ht.PathlineTracker(flow_field)

Next, the starting points of the pathlines need to be defined as a list of arrays, e.g., as follows:

import numpy as np

x0 = np.asarray([0, 5.0, 0])
n = 9
r_arr = np.linspace(0.02,0.2,n)
phi_arr = np.linspace(0, 2*np.pi, n, endpoint=False)

x0_list = []
for r_i in r_arr:
    for phi_i in phi_arr:
        x0_list.append(x0 + np.asarray([r_i*np.sin(phi_i), 0, r_i*np.cos(phi_i)]))

Then the pathlines can be computed using the hemtracer.PathlineTracker.compute_pathlines() method:

pathline_tracker.compute_pathlines(x0_list)

Computing hemolysis

In order to compute the index of hemolysis, we first need to compute the effects of the flow forces on the RBCs. This is done by computing a scalar representative shear rate \(G_s\) along each pathline (see RBC models for details). First, we need to choose a model to compute this shear rate:

cell_model = ht.rbc_model.strain_based.TankTreading()

Depending on the type of model, different configuration options are available, e.g., the sampling rate of the flow gradients for stress-based models and ODE solver settings for strain-based models. For details, refer to the API documentation of the individual models.

In addition, we can choose a model to compute the index of hemolysis from the representative shear rate (see Hemolysis models for details). We will use a power law approach with the correlation from Song et al. [11]:

hemolysis_model = ht.hemolysis_model.PowerLawModel( cell_model, ht.hemolysis_model.IHCorrelation.SONG )

Refer to the documentations of hemtracer.hemolysis_model.PowerLawModel and hemtracer.hemolysis_model.IHCorrelation for details on the configuration options.

Finally, we construct a hemtracer.HemolysisSolver object, which can compute the relevant hemolysis quantities.

hemolysis_solver = ht.HemolysisSolver(pathline_tracker)
hemolysis_solver.compute_representativeShear(cell_model)
hemolysis_solver.compute_hemolysis(hemolysis_model)

Sample Script

The following script shows a complete example of how to use the package:

import hemtracer as ht
import numpy as np
from matplotlib import pyplot as plt
from typing import List, Dict
from os.path import join

def get_seeds():
    """Define seeds for pathlines."""
    x0 = np.asarray([0, 5.0, 0])
    n = 9
    r_arr = np.linspace(0.01,0.2,n)
    phi_arr = np.linspace(0, 2*np.pi, n, endpoint=False)

    x0_list = []
    for r_i in r_arr:
        for phi_i in phi_arr:
            x0_list.append(x0 + np.asarray([r_i*np.sin(phi_i), 0, r_i*np.cos(phi_i)]))
    
    return x0_list

def plot_geffs(hemolysis_solver: ht.HemolysisSolver, models: List[ht.rbc_model.RBCModel], path: str) -> None:
    """Plot Geff for all models in one plot."""
    
    i = 0

    sols = [hemolysis_solver.get_output(model) for model in models]

    for sol_models in zip(*sols):

        plt.figure()
        for sol, model in zip(sol_models, models):
            t = sol['t']
            g = sol['y']
            plt.plot(t, g, label=model.get_name())
        
        plt.legend()
        plt.xlabel("t")
        plt.ylabel("Geff")

        i += 1
        plt.savefig(join(path, "pathline_" + str(i) + "_comp.png"))
        plt.close()

def plot_geff_lagr_eul(hemolysis_solver: ht.HemolysisSolver, model: ht.rbc_model.RBCModel, pathline_tracker: ht.PathlineTracker, path: str) -> None:
    """Plot Geff for Lagrangian and Eulerian data in one plot."""

    i = 0
    geff_lagr_pls = hemolysis_solver.get_output(model)
    geff_eul_pls = pathline_tracker.get_attribute("Geff")


    # Plot pathlines.
    for (geff_lagr_pl, geff_eul_pl) in zip(geff_lagr_pls, geff_eul_pls):

        plt.figure()

        t_lagr = geff_lagr_pl['t']
        g_lagr = geff_lagr_pl['y']
        plt.plot(t_lagr, g_lagr, label="Lagrangian")

        t_eul = geff_eul_pl['t']
        g_eul = geff_eul_pl['y']

        plt.plot(t_eul, g_eul, label="Eulerian")

        plt.legend()
        plt.xlabel("t")
        plt.ylabel("Geff")

        i += 1
        plt.savefig(join(path, "pathline_" + str(i) + "_" + model.get_name() + ".png"))
        plt.close()

def main() -> None:
    # Create flow field object.
    flow_field = ht.EulerianFlowField("flow.vtk")
    flow_field.set_velocity_name("velocity")  # define name of (vector) velocity field, default: "velocity"

    # Load Eulerian data from another file to compare with Lagrangian data.
    flow_field.import_field_data("geff.vtk", field_names_point=["Geff"], field_names_cell=[])

    # Specify which fields are important for the computation.
    flow_field.reduce_to_relevant_fields(["velocity", "Geff", "ref_frame"])

    # Transform data in moving reference frame.
    flow_field.transform_mrf_data(rf_name='ref_frame', rf_rot=2, omega=np.asarray([0, 104.719755, 0]), x0=np.asarray([0,0,0]))

    # Compute pathlines.
    pathline_tracker = ht.PathlineTracker(flow_field)
    pathline_tracker.compute_pathlines(get_seeds(), max_length=100,max_err=1e-5,max_step=0.1,min_step=0.000001,initial_step=0.1,n_steps=10000)

    dt = 0.002          # sampling rate for gradients
    pathline_tracker.interpolate_dv_to_pathlines(sampling_rate=dt) # interpolate velocity gradients to pathlines

    # # Define blood damage models to use.
    model_simp = ht.rbc_model.strain_based.AroraSimplified()
    model_simp.configure_ODE_solver(max_step=dt)
    

    model_tt = ht.rbc_model.strain_based.TankTreading()
    model_tt_corr = ht.rbc_model.strain_based.TankTreadingRotationCorrection()
    model_full = ht.rbc_model.strain_based.AroraFullEig()
    model_full.configure_ODE_solver(max_step=dt)

    model_stress = ht.rbc_model.stress_based.Bludszuweit()
    model_stress.set_sampling_rate(dt)

    correlation_song = ht.hemolysis_model.IHCorrelation.SONG

    cell_models = [model_simp, model_tt, model_tt_corr, model_stress]

    hemolysis_models = [ht.hemolysis_model.PowerLawModel(cell_model, hemolysis_correlation = correlation_song, 
                                    mu = 0.0032, integration_scheme='timeDiff') for cell_model in cell_models]

    hemolysis_solver = ht.HemolysisSolver(pathline_tracker)
    for cell_model, hemolysis_model in zip(cell_models, hemolysis_models):
        hemolysis_solver.compute_representativeShear(cell_model)
        hemolysis_solver.compute_hemolysis(hemolysis_model)
    
    plot_geffs(hemolysis_solver, cell_models, path="../out")


if __name__ == '__main__':
    main()