Brain Deformation From Craniotomy FEM model

[alexrockhill]

I'm attempting to replicate the general procedure in https://doi.org/10.1016/j.media.2019.06.004. The idea is to apply a shift to the boundary condition at a surface. The surface is an ECoG grid that it is apparent is shifted into the brain on a CT compared to the surface of the brain in an MRI before the operation. The voxels of the image can then be shifted appropriately to match this. The code I have so far is below but I'm a bit stuck on implementing this in felupe. Is this possible? And would anyone mind pointing me in the right direction as far as how to do it?

(Requires mne pip install mne pooch nibabel)

import os.path as op
import numpy as np
import matplotlib.pyplot as plt

import mne
import nibabel as nib
import scipy.interpolate

import felupe as fe

misc_path = mne.datasets.misc.data_path()

# get original and projected positions
brainmask = nib.load(op.join(misc_path, 'ecog', 'sample_ecog',
                             'mri', 'brainmask.mgz'))
ras_vox_t = np.linalg.inv(brainmask.header.get_vox2ras_tkr())

raw = mne.io.read_raw(op.join(misc_path, 'ecog', 'sample_ecog_ieeg.fif'))
# two grids, just use one at a time
raw.pick_channels([f'G{i}' for i in range(1, 257)])
trans = mne.coreg.estimate_head_mri_t(
    'sample_ecog', op.join(misc_path, 'ecog'))
montage = raw.get_montage()
montage.apply_trans(trans)  # transform to MRI coordinates

ch_pos_disp = np.array([
    pos for pos in montage.get_positions()['ch_pos'].values()]) * 1000
ch_pos_disp = mne.transforms.apply_trans(ras_vox_t, ch_pos_disp)

raw.info = mne.preprocessing.ieeg.project_sensors_onto_brain(
    raw.info, trans, 'sample_ecog', subjects_dir=op.join(misc_path, 'ecog'))
montage = raw.get_montage()
montage.apply_trans(trans)  # transform to MRI coordinates
ch_pos_orig = np.array([
    pos for pos in montage.get_positions()['ch_pos'].values()]) * 1000
ch_pos_orig = mne.transforms.apply_trans(ras_vox_t, ch_pos_orig)

# fit a surface to the points for deformation
x, y, z = ch_pos_disp.T
y_grid = np.linspace(y.min(), y.max(), 100)
z_grid = np.linspace(z.min(), z.max(), 100)
Y, Z = np.meshgrid(y_grid, z_grid, indexing='xy')

spline = scipy.interpolate.Rbf(
    y, z, x, function='cubic', smooth=5, episilon=5)

X = spline(Y, Z)

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(projection='3d')
ax.axis([0, 256, 0, 256])
ax.set_zlim([0, 256])
ax.plot_surface(X, Y, Z, alpha=0.2)
ax.scatter(*ch_pos_disp.T, c='r')

#########################################################################
# FELupe

aseg = nib.load(op.join(misc_path, 'ecog', 'sample_ecog',
                        'mri', 'aseg.mgz'))
brain = fe.constitution.LinearElastic(E=3000, nu=0.49)
csf = fe.constitution.LinearElastic(E=100, nu=0.49)

# mesh generation
iy, ix, iz = np.where(np.array(brainmask.dataobj) > 0)
nddims = np.array(brainmask.dataobj.shape) + 2

nodemtx = np.zeros(nddims, dtype=np.int32)
nodemtx[ix, iy, iz] += 1
nodemtx[ix + 1, iy, iz] += 1
nodemtx[ix + 1, iy + 1, iz] += 1
nodemtx[ix, iy + 1, iz] += 1
nodemtx[ix, iy, iz + 1] += 1
nodemtx[ix + 1, iy, iz + 1] += 1
nodemtx[ix + 1, iy + 1, iz + 1] += 1
nodemtx[ix, iy + 1, iz + 1] += 1

ndidx = np.where(nodemtx)
points = np.array(ndidx).T
nnod = points.shape[0]

nodeid = -np.ones(nddims, dtype=np.int32)
nodeid[ndidx] = np.arange(nnod)

# generate elements
cells = np.array([nodeid[ix, iy, iz],
                  nodeid[ix + 1, iy, iz],
                  nodeid[ix + 1, iy + 1, iz],
                  nodeid[ix, iy + 1, iz],
                  nodeid[ix, iy, iz + 1],
                  nodeid[ix + 1, iy, iz + 1],
                  nodeid[ix + 1, iy + 1, iz + 1],
                  nodeid[ix, iy + 1, iz + 1]]).T

mesh = fe.mesh.Mesh(points, cells, cell_type='quad')
element = fe.Hexahedron()
quadrature = fe.GaussLegendre(order=1, dim=3)

region = fe.Region(mesh, element, quadrature)
dV = region.dV
V = dV.sum()

displacement = fe.Field(region, dim=3)

u = displacement.values
ui = displacement.interpolate()
dudX = displacement.grad()

F = displacement.extract(grad=True, sym=False, add_identity=True)

kdtree = scipy.spatial.cKDTree(points)

shift = fe.Boundary(displacement, name='shift')
shift.mask = kdtree.query(ch_pos_orig)[1]
shift.points = points[shift.mask]
boundaries = dict(shift=shift)

dof0, dof1 = fe.dof.partition(displacement, boundaries)
ext0 = fe.dof.apply(displacement, boundaries, dof0)

res = fe.newtonrhapson(
    displacement, umat=brain, dof1=dof1, dof0=dof0, ext0=ext0)

[adtzlr]

Hi @alexrockhill,

although the topic is way out of my scope I think I can help you with your FElupe code you provided above. I'll try to explain some things here for you tomorrow.

Best regards,
Andreas

[alexrockhill]

If you have the time, that would be fantastic! I didn't realize that GitHub has a discussion feature instead of issues, thanks for changing it for me to the appropriate label.

[adtzlr]

If you have the time, that would be fantastic! I didn't realize that GitHub has a discussion feature instead of issues, thanks for changing it for me to the appropriate label.

Don't worry, that feature is optional and quite new. 🤠

[adtzlr]

I assume you would like to perform a nonlinear analysis on a mesh with hexahedrons (voxels) in combination with a nearly incompressible material? The following code snippets should do the magic. First I hope you've installed Numba and Pypardiso (if you installed FElupe with all its requirements...) - else, run a terminal and type pip install numba pypardiso. In the following code snippets I use a Hu-Washizu mixed-field approach (u,p,J) to avoid volumetric locking effects. Otherwise, you have to scale down the poisson ratio nu.

In one of my later comments you'll find a link to a shared Google Colab Notebook which contains all your code with my extensions.

[adtzlr]

First, let's try to export the mesh.

mesh = fe.mesh.Mesh(points, cells, cell_type='hexahedron')
mesh.save("mesh.vtk")

[adtzlr]

Both the (though very basic) newtonrhapson function of FElupe and the mixed-field approach only work with hyperelastic materials in FElupe - I suggest you to use the Neo-Hooke material model and convert the Lame constants from E and nu to mu and bulk.

brain = fe.ThreeFieldVariation(fe.NeoHooke(mu=1000, bulk=50000))
csf = fe.ThreeFieldVariation(fe.NeoHooke(mu=100/3, bulk=1650))

[adtzlr]

Next, let's setup the fields. As I already mentioned, we need the cell-wise constant pressure and volume-ratio as additional fields in order to avoid volumetric locking.

displacement = fe.Field(region, dim=3)
pressure = fe.Field(region0)
volumeratio = fe.Field(region0, values=1)

fields = fe.FieldMixed((displacement, pressure, volumeratio))

[adtzlr]

Your approach for the definition of the boundaries was already right - but it may be a bit simplified (and of course, adopted to the mixed-field formulation). In the latest version of FElupe the mask argument is directly available inside the Boundary class. Points are then evaluated automatically by the provided mask. I'm not really sure, do you wish to apply a prescribed displacement on the boundary? Or would you like to create a (stress-like) boundary on the element faces?

Hint: The Boundary class in FElupe only handles prescribed or fixed field values (e.g. displacement components) on mesh points (=nodes).

kdtree = scipy.spatial.cKDTree(points)

shift = fe.Boundary(displacement, mask=kdtree.query(ch_pos_orig)[1])
boundaries = {"shift": shift}

dof0, dof1, offsets = fe.dof.partition(fields, boundaries)
ext0 = fe.dof.apply(displacement, boundaries, dof0)

[adtzlr]

I highly recommend pypardiso as sparse solver because it is magnitudes faster than SuperLU (at least on Windows). The kwargs={"parallel": True} invokes a fast parallel assembly of the stiffness matrix with the help of a JIT-compiled Numba function.

from pypardiso import spsolve

res = fe.newtonrhapson(
    fields, 
    umat=brain, 
    dof1=dof1, 
    dof0=dof0, 
    ext0=ext0, 
    offsets=offsets,
    solver=spsolve,
    kwargs={"parallel": True}
)

fe.save(region, res.x[0], filename="brain.vtk")

[adtzlr]

I've prepared a Google Colab Notebook but unfortunately the RAM is too small. I hope you have more on your local workstation 👍🏻 . Probably the Neo-Hooke material and the mixed field approach is a bit of an overkill - but the code simplifications to a displacement based formulation are straightforward. The LinearElastic material is not very well implemented in FElupe because I use the code mainly for nearly-incompressible hyperelastic simulations on industrial rubber/metal parts.

Good luck, I hope the above code snippets help!

My day is your night 😴
Greetings from Europe/Austria ⛰️

[adtzlr]

Ah, I think I missed something! Do you wish to apply both materials brain and csf to the analysis? If so, better use my other package matADi: pip install matadi. Simply define Neo-Hookean material formulations in MatADi and combine them into a MaterialComposite and finally convert it to a FElupe material.

import matadi
import felupe as fe

brain = matadi.ThreeFieldVariation(
    matadi.MaterialHyperelastic(matadi.models.neo_hooke, C10=500, bulk=50000)
)

csf = matadi.ThreeFieldVariation(
    matadi.MaterialHyperelastic(matadi.models.neo_hooke, C10=100/6, bulk=1650)
)

umat = fe.MatadiMaterial(matadi.MaterialComposite([brain, csf]))

[alexrockhill]

Working through this, it looks like this applies the combined material to all the elements I guess though. If it were possible to add the second material to some of the points, that would make a better model, but maybe first just modeling everything as the same might be easier.

[adtzlr]

That is definitely possible; however, still a bit complicated in FElupe. I'm still trying to find a good codebase for some MultiRegion class. One possible method would be to create two regions with meshes containing the same points but splitted cells. For a mixed-field approach this will get complicated (at least I think so...). The other method is described here.

[alexrockhill]

Wow thanks so much @adtzlr for taking so much time out of your day to help me with this task! I really appreciate it!

To answer your questions:

• The paper describes the displacement as the "loading" of the model for their element-free Galerkin calculation. I'm really not sure the best way to parameterize that into a FEM model, I think with the u vector displacements but my background is much more in signal processing neuroscience so I'm not entirely sure. I can definitely read into it more but, unfortunately, they have not open-sourced their code for the paper I referenced above which makes it a bit trickier, and I think this seems to be the case a lot of this area of biomechanics that the code is not shared. It could be that the element-free Galerkin method works much better since as you mentioned in your first comment, the values of the material properties are not really nearly incompressible but are in fact relatively compressible compared to rubber, steel etc., in which case I may have to just code the algorithm they described from scratch.
• I think starting with the brain material for all of the elements would be reasonable but the aseg file contains labels of the ventricles which would have CSF so elements that fall within that area could be better modeled as csf (4, 43, 14 and 15 are the left and right ventricles, the third and fourth ventricles respectively which contain cerebrospinal fluid (CSF))

I have plenty of RAM on a workstation and I'll execute the code now and post here if I'm able to get it working. Thanks so much again, and good evening from the west coast US :)

[adtzlr]

Hi @recepkaya7! Are you referring to the example above? Could you please describe your problem in a few more words. Sorry, I can't help you with the RAM limit on Google (the code above is quite outdated, FElupe changed a lot in the last year).

[adtzlr]

What do you exactly mean with surface extraction? Obtain points on the surface? Or cells? Or create surface cells?

[recepkaya7]

It's actually like this: I need to apply fem to the brain. When I used the mesh in this example, I had a problem with the newtonrhapson solver due to ram. So I thought of applying fem to a certain part of the brain network. I would be very happy if you could help me with this subject.

[adtzlr]

Okay, I see. But how would you define the boundary condition where your model ends? Treat it as a free surface?

[recepkaya7]

What exactly do you mean when you say from a free surface? And also can't I use boundary conditions in the code ?

[alexrockhill]

This response was very delayed, apologies! (I lost track of quite a few things having a new baby!)

I guess 8 GB of RAM wasn't enough, the process was killed after I tried to run it over night.

I'm not sure where to go from here with felupe, I might try looking into less computationally expensive methods that are less exact like the one in the paper. It seems like a lot to implement, but perhaps I will just give it a go.

Thanks for all your help @adtzlr! You have been more than kind with your time!

[adtzlr]

Thanks, enjoy this wonderful time with as much sleep as possible 😴 ⏰ 🍼.