Inhomogenous anisotropic electrical conductivity

[tcs]

Hi,

I'm attempting to simulate current spread in brain tissue accounting for fractional anisotropy. I have calculated the conductivity tensors from diffusion MRI data, but I'm having trouble working out how to apply them to my model. Anisotropy in the brain varies by region as it is dependent on local white matter tract orientation, so I assume that I will need to input conductivity values that are both inhomogeneous and anisotropic.

Inhomogeneous conductivity data requires input from an analysis output or a S4L cache file. As my conductivity data is not derived from a S4L simulation, does anyone know whether it's possible to import these values?

Thanks!

[PeterStijnman]

Hi,

As far as I know Sim4Life only lets you set the conductivity per voxel, if you want to have anisotropic conductivity you would need to assign different values on the edges of the voxel (on the Yee cell). There might be a possibility to change the dielectric in the hdf file that is created during voxelization step, but my guess is that in there also only the dielectric is saved per voxel and the Yee grid values are only computed once the file is sent to the solver.

So in short you could look into the hdf file that is sent to the solver.

Cheers,
Peter

[tcs]

Thanks for your help. Assigning different conductivity tensors to each voxel is essentially what I am trying to do. Currently, all I can see is that conductivity is assigned to each individual material. Am I looking at this the wrong way?

[PeterStijnman]

In the folder where you have your project saved there should be another folder with the "project_name.smash_results".
In that folder you have potentially 3 files per simulation in the project. It is the UUID of the project and after that the _Voxeler, _Input, and _Output if you ran the simulation. These are .h5 files (hdf). you might find some field in the _Input file that assigns the conductivity to the yee-grid edges (which you can interpret as anisotropic conductivity). If you can change this file to the conductivity you want, and then submit the job from Sim4Life it might work. However, there are then no guarantees that the simulation goes correctly.

Plus the fact that it is highly likely that the conductivity inside this _Input file is still assigned to the voxels themselves and once it is sent to the solver it will compute the conductivity tensor.

If this is the case getting the solver to work with anisotropic conductivity will be very very cumbersome and highly likely will not work correctly (definitely all the things in the Analysis tab will output the correct results, so you would have to export the electric and magnetic fields and do the computations yourself!).

Hopefully, this is a bit of help :)

[montanaro]

Hi yes this is possible in the latest release of Sim4Life, however there isn't a good tutorial for it at the moment (I'll help prepare one). Essentially what is needed is to create a Trivial Producer (an artificial field distribution in the Analysis tab) via Python that has 6 components, representing the symmetric anisotropic tensor values ordered in such a way: xx, yy, zz, xy, yz, xz. Each component should have the full anisotropic values that will be interpolated into the simulation grid. This trivial producer must be assigned the correct units of S/m (this is all handled via a Python script). Once that is done, then in your EM simulation you can assign these values to any given material: Materials -> Electric Conductivity: Anisotropic, Inhomogeneous -> Data Origin: Analysis Output -> Then you can select your field.

I'll hopefully find some time to update this comment here with a sample script.

[montanaro]

@montanaro said in Inhomogenous anisotropic electrical conductivity:

Hi yes this is possible in the latest release of Sim4Life, however there isn't a good tutorial for it at the moment (I'll help prepare one). Essentially what is needed is to create a Trivial Producer (an artificial field distribution in the Analysis tab) via Python that has 6 components, representing the symmetric anisotropic tensor values ordered in such a way: xx, yy, zz, xy, xz, yz. Each component should have the full anisotropic values that will be interpolated into the simulation grid. This trivial producer must be assigned the correct units of S/m (this is all handled via a Python script). Once that is done, then in your EM simulation you can assign these values to any given material: Materials -> Electric Conductivity: Anisotropic, Inhomogeneous -> Data Origin: Analysis Output -> Then you can select your field.

I'll hopefully find some time to update this comment here with a sample script.

Here's a sample script I found to import VTI data and convert it to a Trivial producer in the correct format. Let me know if anything is not clear.

# -*- coding: utf-8 -*-
import os
import numpy as np
import vtk
from vtk.numpy_interface import dataset_adapter as dsa

import XCoreModeling
import XPostProcessor

fp = r"D:\Users\sample\path"
fn = "sample.vti"
f = os.path.join(fp, fn)

reader = vtk.vtkXMLImageDataReader()
reader.SetFileName(f)
reader.Update()

img = dsa.WrapDataObject(reader.GetOutput())
dxx = img.PointData['Dxx']
dyy = img.PointData['Dyy']
dzz = img.PointData['Dzz']
dxy = img.PointData['Dxy']
dyz = img.PointData['Dyz']
dxz = img.PointData['Dxz']

dims = img.VTKObject.GetDimensions()
origin = img.VTKObject.GetOrigin()
spacing = img.VTKObject.GetSpacing()

# scale to meters
scale = XCoreModeling.GetActiveModel().LengthUnits.ToUnit(XPostProcessor.Unit.Meter)
origin = [scale*origin[0], scale*origin[1], scale*origin[2]]
spacing = [scale*spacing[0], scale*spacing[1], scale*spacing[2]]

grid = XPostProcessor.UniformGrid()
grid.SetUnit(0, XPostProcessor.Unit.Meter) # note meters here
grid.SetUnit(1, XPostProcessor.Unit.Meter)
grid.SetUnit(2, XPostProcessor.Unit.Meter)

# create dual grid, so we can assign point data from image as cell centered data
grid.SetDimensions(dims[0]+1, dims[1]+1, dims[2]+1)
grid.SetOrigin(origin[0]-0.5*spacing[0], origin[1]-0.5*spacing[0], origin[2]-0.5*spacing[0])
grid.SetSpacing(spacing[0], spacing[1], spacing[2])
assert grid.NumberOfCells == dxx.shape[0]

sigma = np.zeros([grid.NumberOfCells,6])
sigma[:,0] = dxx
sigma[:,1] = dyy
sigma[:,2] = dzz
sigma[:,3] = dxy
sigma[:,4] = dyz
sigma[:,5] = dxz

aniso_cond = XPostProcessor.DoubleFieldData()
aniso_cond.Grid = grid
aniso_cond.NumberOfSnapshots = 1
aniso_cond.NumberOfComponents = 6
aniso_cond.ValueLocation = XPostProcessor.eValueLocation.kCellCenter
aniso_cond.SetComponentLabel(0, "Dxx")
aniso_cond.SetComponentLabel(1, "Dyy")
aniso_cond.SetComponentLabel(2, "Dzz")
aniso_cond.SetComponentLabel(3, "Dxy")
aniso_cond.SetComponentLabel(4, "Dyz")
aniso_cond.SetComponentLabel(5, "Dxz")
aniso_cond.SetField(0, sigma)
aniso_cond.Quantity.Name = 'Anisotropic Conductivity From DTI (use for Bone)'
aniso_cond.Quantity.Unit = XPostProcessor.Unit('S/m')

producer_t = XPostProcessor.TrivialProducer()
producer_t.SetDataObject(aniso_cond)
producer_t.Description = 'Anisotropic Conductivity'
XPostProcessor.AlgorithmRegistry().AddAlgorithm(producer_t)