From 3D tissue data to impedance using Simpleware ScanFE+IP and COMSOL Multiphysics – a tutorial
Article Category: Tutorial Article
Online veröffentlicht: 26. Apr. 2011
Seitenbereich: 13 - 32
Eingereicht: 05. Apr. 2011
DOI: https://doi.org/10.5617/jeb.173
Schlüsselwörter
© 2011 Fred-Johan Pettersen, Jan Olav Høgetveit, published by Sciendo
This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.
Measurement of electrical bioimpedance is relatively easy to do. It requires little equipment, can be done with small and simple units [1] [2], and due to the simplicity of the equipment, it may be made low-cost. These are all desirable properties in a clinical setting. Electrical impedance in tissues is a function of variables such as ion concentrations, cell geometry, extracellular fluids, intracellular fluids, and organ geometry, so it is easy to imagine that sensitivity to different changes in these variables may be detectable. Unfortunately, it is difficult to establish a link between measured impedance change and a change in one or more of the input variables. This means that electrical bioimpedance measurements, in general, has good sensitivity but poor specificity.
As a result, it may be a challenge to create measurement set-ups of clinical value. For example, if we are to use bioimpedance measurements to diagnose a condition, we need to be sure that our set-up actually detects the given condition. One way to solve this problem is to test a number of set-ups on a number of test-subjects with and without the given condition. However, it is obviously time-consuming for both the persons who do the tests and for the test-subjects, and there may be medical reasons to avoid such tests if not absolutely necessary. And even if we find a set-up that we think is good, we still do not know how much each tissue or organ has contributed to the final measurement result.
One solution to this problem is to create a model where we can see and even quantify the contribution from tissues and organs. This solution is becoming more interesting as computers and software becomes increasingly capable. This tutorial will show how to create such a model and how to perform a virtual experiment on it. The goal is to learn how to model an experiment in order to see what we really measure, and to gain insight in volume impedance measurements in general.
Specifically, this tutorial will take the reader trough the steps involved in making a model and do measurement simulations on it using software tools from Simpleware and COMSOL.
An overview of the tutorial steps are shown in figure 1. It is possible to skip the first parts since deliverables from previous stages are available for download.
Fig. 1
Overview of tutorial. The input data is in DICOM format, and is processed in ScanIP+FE (green). The Model creation and experimenting are done in Multiphysics (blue).
The basis of the model we are going to make is a MRI scan [3, Chapter 12] of a human thigh. We have made the MRI dataset available for download.
The MRI dataset may be seen as a series of greyscale pictures of slices of the thigh. When these 2D pictures are put together, they form a 3D picture where each pixel in a 2D picture is forming a voxel (a volume pixel). The greyscale value of each voxel is saying something about the amount of hydrogen in the tissue; the brighter the voxel is, the more hydrogen there is in the volume the voxel represents. The dataset is stored in a container in a format called DICOM [4] [5]. This makes it possible to use a number of software products to examine the data.
Since we are going to make a model with only a few tissue types, we need to reduce the amount of information in the MRI dataset. The process of replacing the greyscale values with tissue types or organs is called segmentation. Here, we want to segment the data into the most important tissue types: muscle, blood, skin, bone, bone marrow, and fat. Obviously, there are more tissue types present, but we want to create a simple model that will fit into our software tool’s capabilities, so we choose to lump some tissue types together. After segmentation we have a 3D dataset where each voxel has a name or number that is representing a tissue type. Figure 2 shows a typical SIP window. When the segmented model is made, we can create copies of it and add electrodes.
Fig. 2
A typical SIP window. In this particular case, after segmentation.
It is also possible to use datasets that are pre-segmented. An example of a pre-segmented dataset is the VOXEL-MAN dataset [6]. In such a case, the pre-segmented dataset must be imported and treated just like the dataset we segment ourself.
Our simulation tool is not capable of using a segmented model directly, so we have to create a mesh. A mesh is another representation of the data where the geometric information is represented by points in 3D space, and lines between these points. Figure 3 shows an example geometric model seen from the outside. Figure 4 shows the same with the interior features made visible. It is possible to represent this with a mesh as shown in figure 5. The mesh may be coarse (figure 5) or fine (figure 6). The simulation tool will then use the mesh to calculate physical state variables (scalars) for the the intersection points in the mesh and vectors along the lines between the points. This method is also known as Finite Element Method (FEM) [7]. If a fine mesh is used the simulation tool has more data to do calculations on, and simulations will therefore take longer time and require more computer memory.
Fig. 3
Example geometric model.
Fig. 4
Example geometric model with interior features shown.
Fig. 5
Mesh representation of example geometric model. Coarse mesh.
Fig. 6
Mesh representation of example geometric model. Fine mesh.
When a mesh is created, we import it into our simulation tool and merge it with electrical properties of the tissues and materials used to create a model that is representing both geometrical and electrical properties of the thigh. The electrical properties we use are frequency dependent electrical conductivity (
Fig. 7
Electrical properties for the tissues used in this tutorial [8]. The properties are electrical conductivity (
Fig. 8
A typical MPH window with a simulation result in the graphics-pane.
To be able to perform all steps, these resources are required: Scan IP+FE (SIP) from Simpleware [9] , version 4.2, Multiphysics, version 4.1 (MPH) from COMSOL [10] with the AC/DC Module. A computer capable of running the above software. A minimum of 8 GB RAM is recommended.
SIP is used for segmentation and meshing, while MPH is used for experiment set-up and simulation.
The files generated by SIP are made available, so it is possible to perform the MPH steps even if you do not have SIP.
A number of files are made available for this tutorial. These include raw data, mesh files, tissue data, and MPH model files in various degree of completion.
The supplied data files for SIP are:
DICOM image stack of a two thighs. Input to SIP. DICOM is the name of a data-wrapper for radiological images. Many CT [11] and MRI [12] systems produce DICOM datasets [4].
Presegmented data of a thigh.
SIP project with electrodes.
Spreadsheet file with statistics from SIP. Column with MPH-domain mapping is added.
The supplied data files for MPH are:
MPH mesh generated by SIP (form tutorial.sip).
Text file containing expressions for
MPH file with materials pre-defined.
MPH file with materials and mesh. Experiment set-up is not defined.
MPH file with materials, mesh, and experiment set-up.
MPH file with finished tutorial. This file contain everything including results.
The files are packaged in a ZIP archive made available at:
If you are not familiar with bioimpedance and it’s principles or the software used in this tutorial, the following reading material is recommended:
Bioimpedance in general:
ScanIP:
COMSOL Multiphysics:
An overview of the process of creating a mesh form a MRI dataset is shown in figure 9.
Fig. 9
Overview over the segmentation and meshing meshing process. The raw data is coming from a DICOM source such as MRI or CT. As an alternative, one may use pre-segmented data (dashed arrow). This is done using Scan IP+FE, and the output is a mesh usable by Multiphysics.
If your raw-data is coming from a MRI machine or a CT, the first step is to import and segment the data. We will use MRI-images of a thigh in this tutorial.
Import of the supplied DICOM dataset is described below:
Start ScanIP.
Select File
Browse to the directory MR femur. Do not go into the underlying folders.
Select the uppermost series, and click OK.
In Window and Level, click in the histogram charting area and then select Full Data Range.
In Crop, play around with the slides together with the Preview-slide. The Crop settings used here are:
Lower X: 24, higher X: 439.
Lower Y: 235, higher Y: 583.
Lower Z: 0, higher Z: 32.
Click OK. There will be a warning about non constant Z spacing. This will be ignored in this tutorial.
7. Save the file.
Use the sliders on the three 2D windows in ScanIP to examine the dataset. Notice that since this is MR data, the bones are black. Slice 12 should look like figure 10.
Fig. 10
Slice 12 of imported MR data.
The bone (Femur) is the black ring in the middle of the thigh. We are going to use a tool called
Display slice 12, (
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Number of iterations: 3.
Multiplier: 2.5.
Initial neighborhood radius (pixels): 1.
Mode: 3D.
Create new mask.
Perform on all slices.
Click at position (154, 170, 12). Se pointer position in the lower left corner of the ScanIP-window.
Rename the created mask. Right-click the mask and select Rename.... Enter the text Bone.
Change colour. Right-click the mask and select Change colour
In the 3D-window select Fast Preview and click Refresh. This should show a bone that looks very much like a femur. If not, choose Undo and try again.
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 3.
Y radius: 3.
Z radius: 1.
Sub-voxel component (pixels): 0.000.
Filter: Close.
Click Apply. The result should look like figure 11.
In the 3D-vindow click Refresh, and you should get something like figure 12. You have to manipulate the 3D-view using he mouse to get an identical view.
Fig. 11
Slice 12 of segmented bone.
Fig. 12
3D-view of segmented bone.
The bone marrow is the gray tissue in the middle of the bone. We are going to use a tool called
Display slice 12, (
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Number of iterations: 3.
Multiplier: 3.00.
Initial neighborhood radius (pixels): 1.
Mode: 3D.
Create new mask.
Perform on all slices.
Click at position (170, 150, 12).
In the 3D-window select Fast Preview and click Refresh. This should show a bone with bone marrow. If not, choose Undo and try again.
Rename the mask to Bone marrow, and change colour to yellow.
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 3.
Y radius: 3.
Z radius: 1.
Sub-voxel component (pixels): 0.000.
Filter: Close.
Click Apply.
Select Current tool
Apply on active mask.
X radius: 1.
Y radius: 1.
Z radius: 1.
Sub-voxel component (pixels): 0.000.
Filter: Dilate.
Click Apply.
In the 3D-window click Refresh.
We need to address the gap between the bone marrow and the bone. The chosen method is to fill the bone (make it non-hollow), and place the bone marrow mask above the bone mask. This way, the bone marrow will overwrite the bone in the meshing process, and the bone will be meshed as a hollow bone with bone marrow inside it.
Select Bone mask.
Display slice 12, (
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 1.
Y radius: 1.
Z radius: 1.
Sub-voxel component (pixels): 0.000.
Filter: Dilate.
Click Apply.
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Apply from active mask.
Mode: 3D.
Merge with mask.
Perform on all slices.
Click at position (170, 150, 12).
Select Current tool
Apply on active mask.
X radius: 1.
Y radius: 1.
Z radius: 1.
Sub-voxel component (pixels): 0.000.
Filter: Erode.
Click Apply.The result should look like figure 13.
If the bone marrow mask is below the bone mask, drag the bone marrow mask onto the bone mask (it will now appear above the bone mask in the list).
Fig. 13
Slice 12 of segmented bone and bone marrow. Note that it may be difficult to see that the colours of bone and bone marrow are different in grayscale prints.
The
The blood vessels are small on the MR dataset, from 1 to around 8 pixels wide. And the since the voxels are 7 mm long in the Z-direction, it is practically impossible to model blood vessels that do not go more or less parallel to the Z-direction since the voxels are not close to cubic in shape. We are going to use a tool called
Display slice 12, (
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Number of iterations: 2.
Multiplier: 1.5.
Initial neighborhood radius (pixels): 1.
Mode: 3D.
Create new mask.
Perform on all slices.
Click in the middle of the dark blood vessel at location (337, 192, 12).
Rename the new mask to Blood, and change the color to Artery.
Change the Mask operation to Merge with mask.
You may add more blood vessels by clicking on these locations: (354, 190, 12) and (342, 138, 15) (note that the last coordinate requires change of slice to 15), or you may play around with the settings and where you click. If your segmentation should produce really bad results you can always use undo.
In the 3D-window click Refresh to see something similar to figure 14.
Fig. 14
3D-view of bone, bone marrow, and blood.
It is smart to use the Refresh-button often in the 3D-view.
The fat and connective tissue is the bright parts on the outside of the muscles and bone. For simplicity both fat and connective tissue is grouped into fat. We are going to use a tool called
Display slice 12, (
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Number of iterations: 3.
Multiplier: 3.
Initial neighborhood radius (pixels): 1.
Mode: 3D.
Create new mask.
Perform on all slices.
Click in the middle of the bright area at location (111, 262, 12).
Rename the new mask to Fat, and change the color to Fat.
In the 3D-window click Refresh. You will see that the mask is not completely defining the tissue, so we need to close some holes.
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 3.
Y radius: 3.
Z radius: 1.
Sub-voxel component (pixels): 0.000.
Filter: Close.
Click Apply.
In the 3D-window click Refresh. There are still some open holes that can be filled semi-manually:
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Apply from active mask.
Mode: 2D
Merge with mask
Perform on active slice.
Go through all slices and click in uncovered areas that should be marked as fat. Use undo if the algorithm is marking too much. It is OK that this mask is overlapping other masks as it will be given least priority before meshing.
The muscle tissue is basically what is left. This makes it easy to segment it using boolean techniques.
Right-click on the Fat mask, and select Duplicate.
Rename the new mask to Muscle and give it the color Muscle.
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Brush shape Disk.
Brush size 3.
Smooth edges.
Paint operation: Paint.
Perform on active slice.
On slices 0 through 7 paint the gap around (270, 305, x).
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Apply from active mask.
Mode: 3D
Merge with mask
Perform on all slices.
Click inside a muscle area to make the muscle mask cover about everything in the thigh.
Right-click the Muscle mask and select Boolean operations
Right-click the Muscle mask and select Boolean operations
Right-click the Muscle mask and select Boolean operations
Right-click the Muscle mask and select Boolean operations
13In the 3D-window click Refresh
The skin tissue is a layer that we decide is 1 mm thick (i.e approximately 2 pixels), and is covering the outside of the thigh. This is also a good candidate for boolean operations.
Right-click the Fat mask, and select Duplicate.
Rename the new mask to Skin and give it an off-white colour.
Right-click the Skin mask and select Boolean operations
Right-click the Skin mask and select Boolean operations
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 2.
Y radius: 2.
Z radius: 1
Sub-voxel component (pixels): 0.000.
Filter: Dilate.
Click Apply.
Note that the Skin-mask is actually covering everything in to the bone. This is done deliberately since we need it for a later boolean operation. It will not affect our model since the precedence of the Skin mask is the least of all masks and the other masks will overwrite the Skin mask where appropriate.
It is necessary to expand the Fat mask so that it includes Fat, Muscle, Blood, and Bone. This is to make sure there are no holes in the structure. This is particularly important when we apply filters on the dataset. Since there is a mask priority system in SIP these operation will not cause the Fat mask to overwrite other masks.
Right-click the Fat mask and select Boolean operations
Right-click the Fat mask and select Boolean operations
Right-click the Fat mask and select Boolean operations
Re-arrange the mask-list so that this is the new order from the top: Bone marrow, Bone, Blood, Muscle, Fat, and Skin.
Save your work.
The segmentation is now finished, and the result is saved in a SIP project file. Se figures 15 and 16 for how it should look like. This file will now be used for meshing. It is recommended to make a back-up copy of it.
Fig. 15
Slice 12 after segmentation.
Fig. 16
3D-view of thigh after segmentation.
In some cases we have pre-segmented datasets available, and we may choose to work on these. An example of such dataset is the VOXEL-MAN [6]. The VOXEL-MAN dataset is not directly importable, but with a few modifications such as merging organs where these are divided, it can be used. Import of pre-segmented data is straight-forward in SIP by using File
In principle, it is possible to use the segmented dataset directly to create a mesh, but applying some filters and resampling might improve the quality of the mesh and reduce the size of the resulting mesh.
Before filtering is started, it is important to understand that the order of masks is of importance. The masks may overlap each other, partly or completely, and if so happens, the upper mask (in the dataset browser in SIP) is effectively overwriting the lower masks. So if a muscle mask and blood vessel mask is overlapping, the blood vessel will be overwritten by the muscle mask if the muscle mask is above the blood vessel mask. To rectify this, simply drag the blood vessel mask over the muscle mask. Although overlapping masks might be a problem, the mask importance can be useful for us. Consider the case where there are several small holes in our mask set caused by filtering. Now we can simply define a mask that covers all those holes and call it
The first task is to resample the dataset. This may not be strictly necessary given enough computer memory and time to simulate, but we want to have a smaller dataset and we want to demonstrate down and up-sampling. The cost is poorer resolution.
Select Data
Settings to use:
Units: Percent change (%).
X: 30.
Y: 30.
Z: 150.
Background Interpolation method: Majority wins.
Mask Interpolation method: Majority wins.
These settings will reduce the number of voxels in X and Y direction, and increase the number of voxels in Z direction.
Click OK.
The number of voxels are reduced to approximately 13 % to around 620,000 voxels. This causes loss of data, but makes it possible to create a relatively small model that will fit in memory and simulate in a reasonable time.
If we look in the 3D-window we see that the skin mask has holes. These holes are difficult to fix using filters, so we are going to re-create the skin mask. Keep in mind that the skin is difficult to model since it is a very thin structure. There are two requirements to the skin: It must be possible to smooth it using filters without creating new holes, and it must be thick enough (in terms of voxels) to create a good quality mesh. A skin that is 3 mm (2 voxels) thick will fulfill our requirements. Having a skin that might be thicker or thinner than the real skin is acceptable for this particular case since the skin does not contribute much to a four-electrode measurement.
Right-click the Skin mask and select Delete.
Right-click on the Fat mask, and select Duplicate.
Rename the new mask to Skin and give it an off-white colour.
Drag masks so that the skin mask is at the bottom.
Right-click the Skin mask and select Boolean operations
Right-click the Skin mask and select Boolean operations
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 2.
Y radius: 2.
Z radius: 1.
Sub-voxel component (pixels): 0.000.
Filter: Dilate.
Click Apply.
The filtering we are going to do is just basic smoothing on all masks. It is recommended to hide the masks (make them invisible) you are not working on and to Refresh the 3D-view after each operation.
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 3.
Y radius: 3.
Z radius: 3.
Sub-voxel component (pixels): 0.000.
Apply the filter with these settings to these masks: Bone marrow, Bone, Muscle, Fat, and Skin.
Select Current tool
Apply on active mask.
X radius: 4.
Y radius: 4.
Z radius: 4.
Sub-voxel component (pixels): 0.000.
Filter: Close.
Make Blood active and click Apply. Note that the blood vessels are not preserved perfectly. This is due to the small feature size. It is possible to improve the blood vessels by manually editing the masks and by further experimenting with the filters. If deemed interesting, the blood vessel modeling should be done more thoroughly.
Save your work.
Figure 17 shows how the filtered masks are looking in a 3D-view. Use the visualization options to change opacity of masks.
Fig. 17
3D-view of thigh after filtering.
More on resampling
It is easy to loose vital information when resampling since we, in this case, are reducing the number of voxels. It may be necessary to handle some masks in a special way when this is the case. If we consider the case where a relatively thin feature like a blood vessel should be preserved, a resampling will probably fragment this. One way to handle this is:
Decide target voxel size. It should be small enough to represent the smallest feature properly.
Resample to the largest voxel size where your small feature still is conserved. This is probably somewhat smaller than smallest target voxel size.
Save a copy of your project in a temporary file.
Continuing in your main project file, resample to the target voxel size.
Resample to the voxel size of your temporary file.
Import the mask(s) from your temporary file. Your current masks will not be overwritten.
Make the mask representing the small feature active.
Right-click the mask and select Boolean operations
The voxel size is now smaller than the target voxel size, but this should not be a problem since the feature size is the same as the target voxel size. The meshing tool will still be capable of creating a mesh that is small. This method may require some padding of the dataset to make it possible to resample and merge two datasets. Filtering of the small features may be needed.
The dataset is now segmented and filtered, and is ready for the final step before meshing which is adding electrodes. It may be a good idea to save a copy of the model now so it is possible to go back and try other electrode configurations if wanted. We will add four electrodes here. The general idea is to draw electrodes using the segmentation tool Paint. The Electrodes will be drawn semi-invasive since this is easy to do, and then the invasive part will be removed using boolean operations.
Create electrode mask: Right-click over the mask list and select Create new mask.
If necessary, move the mask to the top of the list.
If the mask is invisible, toggle visibility.
Change name to Electrodes and colour to Blue.
In the Tool selector, click on the Segmentation-tab.
Select Current tool
Brush shape Square.
Brush size 9.
Smooth edges.
Paint operation: Paint.
Perform on active slice.
Draw electrodes. Z-value for the electrodes are [7-9], [15-17], [30-32], and [38-40], X-value is 67. Place the cross in the paint cursor just on the outside of the skin and click to paint.
Make sure the electrode mask is active.
In the Tool selector, click on the Basic filters-tab.
Select Current tool
Apply on active mask.
X radius: 3.
Y radius: 3.
Z radius: 3.
Sub-voxel component (pixels): 0.000.
Click Apply.
Right-click the mask and select Boolean operations
Figure 18 shows how the thigh with electrodes in a 3D-view.
Fig. 18
3D-view of thigh with electrodes.
Meshing
The dataset is now ready for meshing. The tasks involved in the meshing are affected by the type of tool and type of simulation we plan to do. Here, our target is finite element (FE) simulation using COMSOL Multiphysics. Here is how to create a mesh:
Create Model configuration: Right-click on Masks and select Create new FE model.
Select all masks (click on the top-one and shift-click on the last).
Right-click the masks and select Add to active model.
Right-click the model (Model 1) and select Model configuration....
Enter these settings:
Model type: FE.
Export type: Comsol volume.
Pre-processing
Pre-processing
Pre-processing
Pre-processing
Volume meshing
Volume meshing
Click Close.
In the 3D view, select Models and Fast preview. Click Refresh. The thigh should not have changed.
In the 3D view, select Full model, and click Mesh. This operation may take some time.
In the 3D view, click Export, and enter a file name for the mesh.
Click OK and Close in the relevant boxes after export is done.
If the Project information toolbox is not open, click Toolboxes
Click Update Stats and Export Data, and enter a file name.
Save the ScanIP file.
Open the file generated in previous step in a spreadsheet program such as OpenOffice Calc or similar.
Add a column between the two first columns, and give it the name MPH domain. Give each tissue a number. The tissue at the bottom of the mask list in SIP should have 1, the tissue above should have 2, and so on. This is necessary since the exported data is not sorted by mask precedence as we would like it to be.
Sort the tissue rows by ascending numbers in the new column if needed, and save the file for later use.
The exported mesh is now ready to be imported into COMSOL Multiphysics (MPH).
Models are created using two inputs: Tissue data and mesh. See figure 19 for an overview of the process.
Fig. 19
Overview of modeling process. The inputs are tissue data and the previously generated mesh, which are being merged to create a model. This is the base for our bioimpedance experiment. This is done using Multiphysics.
MPH is material-oriented. That means that each model must have one or more materials with the relevant material properties defined. We only need two parameters: Electrical conductivity (
Since the tissue parameter library is ready-made, the next step is to create the materials. Here is how to create a model file containing only the materials we need:
Start MPH.
In the Model Wizard, select 3D and click finish button. We are not using the Model wizard to add physics and studies sine it necessary to learn to do this.
In the Model Builder, right-click the Materials-node, and select Material to add a new material.
Right-click the new material and rename it to SkinDry.
Click the new SkinDry node.
In Settings, expand Material
Right-click Electric Conductivity, and select Add to Material.
Right-click Relative Permittivity, and select Add to Material.
In the value-fields below, enter these expressions for sigma and epsilonr (keep in mind that MPH is case-sensitive):
SkinDry sigma
SkinDry epsilonr
Repeat steps 3-10 for FatAverageInfiltrated, Muscle, Blood, BoneCortical, and BoneMarrowInfiltrated.
The materials we need are now created, but we need to load the expressions for the conductivities and permittivities from the file in the supplementary data. Here is how to import the tissue parameter library:
In the Model Builder, right-click Global Definitions and select Parameters.
In the Settings
Locate and select the file tissue parameter library.txt.
This action will add all tissue data available from [8]. Many of these will not be used in this tutorial, so you may remove the unused parts from this text-file if you prefer. We also need to create a metal for our electrodes. For this, we will use a built-in material:
Select the Material Browser-tab.
Expand the Built-In-node.
Scroll down to and right-click Steel AISI 4340.
Select Add Material to Model.
This material has all necessary parameters defined, so we do not have to do anything more with it.
A MPH file containing the materials is available in the supplementary data. It is named tutorialmaterials.mph.
The mesh created using ScanIP needs to be imported. This is the procedure:
Right-click the Mesh-node and select Import.
In the Settings-pane, select Browse and find your mesh-file. If you prefer to use the supplied mesh, it is called tutorial.mphtxt.
Only Split on material data should be checked.
Click Import. This operation may take some time.
You should see the mesh in the Graphis-pane. You will have to zoom in to see the mesh details.
The final step in creating a model is to connect mesh data and material data. Since MPH is material-orientated, we go through our materials and add domains to it. If transparency is not turned on in the Graphics pane, then click on the Transparency button. Here is the procedure:
Open the spreadsheet created when exporting the mesh in OpenOffice Calc or similar. If you are using the supplied files, this is called tutorialStats.csv.
In MPH, expand the Materials-node.
Select SkinDry, the first material in the list.
In the Settings-pane, change Selection to Manual.
Select all domains and click Remove from Selection (the minus-button).
If the Selection List-tab is not available (it is normally placed next to the Model Builder-tab), Select Options
In the Selection List, select the number corresponding to SkinDry according to the spreadsheet. It is 1 for the supplied data. This should turn the skin in the model red.
Click Add to Selection (the plus-button). This should turn the skin in the model blue.
Select the Model Builder-tab.
Select the next material, FatAverageInfiltrated. From this material an onwards, the selection is already set to manual, so it is no need to change that.
Select the Selection List-tab.
In the Selection List, select the number corresponding to Fat according to the spreadsheet. It is 2 for the supplied data.
Click Add to Selection (the plus-button).
Repeat 9 - 13 for the rest of the materials. The connections (for the supplied data) are shown in the table 1.
Tissue mapping table for mapping data between the ScanIP mesh and MPH.
| ScanIP name | MPH name | MPH Domain |
|---|---|---|
| Skin | SkinDry | 1 |
| Fat | FatAverageInfiltrated | 2 |
| Muscle | Muscle | 3 |
| Blood | Blood | 4 |
| Bone | BoneCortical | 5 |
| Bone marrow | BoneMarrowInfiltrated | 6 |
| Electrodes | Steel AISI 4340 | 7 |
Figure 20 shows the model ready for experimenting.
Fig. 20
3D view of the model with electrodes selected.
It is now time to save your file. A model file containing the model is available in the supplementary data. It is named tutorial-model.mph.
There are a number of experiments that can be done. The one presented in this tutorial is a basic four electrode transfer impedance measurement. See figure 21 for an overview.
Fig. 21
Overview of experiment processes. The process is based on the model previously defined. This is done using Multiphysics.
In this experiment, we will use our ready-made model to measure transfer impedance in a four-electrode system like the one shown in figure 22. The measured transfer impedance is what we can reproduce using suitable equipment. We will also have a look at some parameters that we cannot measure in a simple way. These are: Current paths, sensitivity, volume impedance density. Furthermore, we will calculate how much the different tissues contribute to the measured transfer impedance. These last parameters require a measurement set-up like the one in figure 23 where we explore the principle of reciprocity.
Fig. 22
Normal 4-electrode transfer impedance measurement set-up with current stimulus on the outer electrodes and voltage measurement on the inner electrodes. This is a real-world set-up.
Fig. 23
Measurement set-up for exploring the principle of reciprocity. This set-up has current stimuli on both outer and inner electrode pairs. This set-up is not used for direct measurements.
First of all we need a model to experiment on. If you have not made one so far in the tutorial, make a copy of tutorial-model.mph and give it a new name of your choice. Open this file in MPH.
First, add two physics interfaces. The first of these (
Right-click the Model-node, and select Add Physics.
In the Model Wizard-pane, expand AC/DC.
Select AC/DC
Click Add Selected (the plus-button) twice.
Click the Finish-button (the flag).
Since we are working with a linear model we may just as well use 1 A currents. We are specifying the currents in the physics interfaces by adding a Terminal and a Ground-node in each physics interface. Here is how to do that:
Right-click the first physics interface (
In the Settings pane set the terminal type to Current, and set
Add the positive CC electrode (CCP on figure 22) to the terminal using one of two methods:
Click on the Selection List-tab.
Try clicking on the numbers until the top CC electrode turns red (the electrode closest to the thick end of the thigh). In the supplied model, it is number 11.
Click Add to Selection.
Or:
Click on the top electrode in the Graphics-pane until it turns red.
Right-click anywhere in the Graphics-pane. Note that the selection is turning blue, and that a number appear in the Boundaries list.
On this particular electrode it is possible to select the sides of the electrode (number 317) as well, but it is not necessary since the conductivity is very high in the electrode material.
Right-click the first physics interface (
Add the other CC (CCN on figure 22) node to Ground using one of two methods:
Add the negative CC electrode to the terminal using one of two methods:
Click on the Selection List-tab.
Try clicking on the numbers until the bottom CC electrode turns red (the electrode closest to the thin end of the thigh). In the supplied model, it is number 24.
Click Add to Selection.
Or:
Click on the bottom electrode in the Graphics-pane until it turns red.
Right-click anywhere in the Graphics-pane.
Repeat the above steps for the second physics interface (
Next, we must add a study so we can tell MPH what to do.
Rigth-click the top node (
In Select Study Type, select Frequency Domain and click the Finish-button.
In Settings-pane
In Physics Selection, make sure that only
This is enough information for MPH to do a simulation. Out of this simulation we can get transfer impedance, current densities, potentials and so on. But we can not get sensitivity, volume impedance density, or information on how much a given part is contributing. To get this we must set up another study step:
Rigth-click the Study-node and select Study Steps
In Settings-pane
In Physics Selection, make sure that only
Since the initial values for the solver is 0 V for all points in the model, and the solution may be far from 0 V, there is a fair chance that the solver will have some problems. The solver may not be able to find a solution, the solution may be wrong, or the solver may use long time to reach a solution.
One method of helping the solver is to specify Initial Values for
Fig. 24
Initial voltage values for simulation. The z-axis is going along the thigh, while the y-axis is going from front to back through the thigh.
Since voltage is the dependent variable when using the
Expand the first physics interface node (
In the Settings-pane, Initial Values enter this expression:
750+(750*exp(-(((z-0.18)ˆ2)/(2*0.01ˆ2)))-750*exp(-(((z-0.03)ˆ2)/(2*0.01ˆ2))))*(1/(1+100*y))
Expand the second physics interface node (
In the Settings-pane, Initial Values enter this expression:
750+(750*exp(-(((z-0.15)ˆ2)/(2*0.01ˆ2)))-750*exp(-(((z-0.06)ˆ2)/(2*0.01ˆ2))))*(1/(1+100*y))
Solving, or running a simulation, is done by:
Select the Study-node.
Click the Compute-button (the equal-sign).
Note that simulations may take some time and that the memory requirements is close to 4 GB. It may be necessary to close all other programs running on the computer to free up memory. It may also be a good idea to clean up MPH memory usage by saving your model, close and re-start MPH, and re-load your model before solving.
In some cases, the mesh generated by ScanIP+FE has elements with low quality. There are several reasons for this, but generally spoken, increased mesh-quality requires either more work and/or larger mesh (in other words: finer mesh, see figures 5 and 6). If the mesh has low quality, there will be warnings about this in the solver nodes in MPH (Study
It is possible to check the quality of the mesh by right-clicking the Mesh-node and select Statistics. The most important number is Minimum element quality, which should be above 0.1. The mesh used in this tutorial is not as good as desired, but the amount of low-quality mesh-elements is relatively low, so it is accepted. To see where the mesh has low-quality it is possible to plot the low-quality elements of the mesh. Here is a description of how:
Right-click the Results node, and select 3D Plot Group.
Right-click the 3D Plot Group, and select Mesh.
In the Settings-pane expand Element Filter.
Check Enable filter.
Enter the expression
Click the plot-icon at the top right of the Settings pane.
This procedure will display all parts of the mesh where the quality is less than 0.5. By changing the number in the expression to 0.1, it becomes clear that the overall quality of the mesh is rather good even if some elements are poor.
If you choose to work with a mesh containing low-quality elements, it is possible to look away from these nodes in calculations of sensitivity, volume impedance density and derived quantities. This is done by adding a
The results of a simulation is placed in the Data Sets-node (which is placed in the Results node). Each study gets a solution node here. For this experiment, there are only one node since we have only one study with one study step. There are a number of parameters that can be visualized in the results part, but here, only a few relevant will be demonstrated.
MPH will create a few 3D Plot Groups, but we will delete these and create our own to learn how to do it:
Expand the Results node, and right-click the first 3D Plot Group.
Select Delete, and confirm that you want to delete the plot group.
Repeat for the remaining 3D Plot Groups
The first we want to see is surface potential resulting from the current injected in the CC electrodes. Here is how to create this view:
Right-click the Results node, and select 3D Plot Group.
Right-click the 3D Plot Group, and select Rename. Give the plot group a name such as ec surface potential.
Right-click the 3D Plot Group, and select Surface.
In the Expression field, enter V. V is the free variable for
Click the plot-icon at the top right of the Settings pane. This should generate a plot as shown in figure 25.
Fig. 25
Surface plot of potential generated by injecting 1 A through the CC electrodes.
Next we want to see potential inside the model. To do this we use a slice plot. Here is how to create this view:
Right-click the Results node, and select 3D Plot Group.
Right-click the 3D Plot Group, and select Rename. Give the plot group a name such as ec slice potential.
Right-click the 3D Plot Group, and select Slice.
In the Expression field, enter V.
In the Plane Data section, enter
Plane type: Quick
Plane: yz-planes
Entry method: Coordinates
x-coordinate: 0.103
Click the plot-icon at the top right of the Settings pane. This should generate a plot as shown in figure 26. Feel free to play around with the other settings, for example change the slice planes.
Fig. 26
Slice plot of potential generated by injecting 1 A through the CC electrodes.
The surface and slice plots are useful for potentials, but if we want to have a look at vector quantities, we need to use other plots such as the arrow volume plot. The next quantities we are going to investigate are current densities for both the normal current that is applied in the CC electrodes and the reciprocal currents that is applied in the PU electrodes. We are going to investigate both in the same plot, but with different colours. Here is how to create this view:
Right-click the Results node, and select 3D Plot Group.
Right-click the 3D Plot Group, and select Rename. Give the plot group a name such as current densities.
Right-click the 3D Plot Group, and select Arrow Volume.
The Expression fields may be filled out, but in case it is not, click the Replace Expression icon and select Electric Currents
Change Coloring and Style
Change Coloring and Style
Click the plot-icon at the top rigth of the Settings pane to see the normal current densities.
Right-click the 3D Plot Group, and select Arrow Volume.
In the Expression field, click the Replace Expression icon and select Electric Currents 2
Change Coloring and Style
Change Coloring and Style
Click the plot-icon at the top rigth of the Settings pane to see both the normal and reciprocal current densities.
Click the Zoom In button once.
This should generate a plot as shown in figure 27. Depending on the default values, it may be a good idea to play with the Arrow Positioning and Coloring and Style. Note that it is possible to use slice plot (as described earlier) to plot the absolute value (no direction information) of the current densities. The way to do this is to create a slice plot and use the expression ec.normJ.
Fig. 27
Plot of current densities for normal (red) and reciprocal (blue) currents.
To examine sensitivity and volume impedance density we must have a look at the expressions for the principle of reciprocity [13], which are
Or as an expression in MPH:
(ec.Jx*ec2.Jx+ec.Jy*ec2.Jy+ec.Jz*ec2.Jz) /((1[A])ˆ2). The expression for volume impedance density is
Or as an expression in MPH: (ec.normE/ec.normJ)*
((ec.Jx*ec2.Jx+ec.Jy*ec2.Jy+ec.Jz*ec2.Jz) /((1[A])ˆ2)).
It is important to keep in mind that the sensitivity (S) is not telling anything about what we measure in our physical measurement set-up, it merely tells us how much weight the local (voxel) impedance is getting. We have to look at volume impedance density (VID) to find the local contribution to the measured transfer impedance.
Since S and VID are scalars, we choose to use slice plots to visualize them.
Right-click the Results node, and select 3D Plot Group.
Right-click the 3D Plot Group, and select Rename. Give the plot group a name such as sensitivity.
Rigtht-click the 3D Plot Group, and select Slice.
In the Expression field, enter the MPH-version of equation 2.
Change Plane Data
Click the plot-icon at the top rigth of the Settings pane. The result should be something like figure 28.
Right-click the Results node, and select 3D Plot Group.
Right-click the 3D Plot Group, and select Rename. Give the plot group a name such as volume impedance density.
Right-click the 3D Plot Group, and select Slice.
In the Expression field, enter the MPH-version of equation 3.
Change Plane Data
Click the plot-icon at the top rigth of the Settings pane. The result should be something like figure 29.
Fig. 28
Slice plot of sensitivity inside the thigh.
Fig. 29
Slice plot of volume impedance density inside the thigh.
Feel free to play around with the other settings, for example change the slice planes.
Figure 29 clearly show that there are a lot of things happening around the electrodes, and that there are large negative and positive contributions. If we zoom in on the electrodes wee will clearly see the layers of tissue.
We are now going to look at some results that are represented by scalar quantities instead plots. Let us have a look at how much the bone, and how much the muscle tissue contribute to the measured results. This is particularly interesting if we want to target one particular tissue or organ. This is how:
Define an integration operator for the complete model and for the tissues or organs in question.
Expand your model node.
Rigth-click Definitions, and select Model Couplings
Add the muscle domain. Select it either by clicking at the model in the Graphics pane, or by using the Selection list pane. If you use the supplied files, it is domain number 3.
Repeat for the bone domain. If you use the supplied files, it is domain number 5.
Repeat for all domains. By all domains we mean that all domains should be selected. This is easiest done by changing Selection from Manual to All domains.
Repeat for the skin domain. If you use the supplied files, it is domain number 1.
You should have four integration operators now. The next step is to add expressions for calculation of the desired quantities.
Expand your model node.
Right-click Definitions, and select Variables.
In the variables list add four variables named
Fill out the expression-field using
intop1(VID)
intop2(VID)
intop3(VID)
intop4(VID)
where VID is the expression in equation 3. It is possible to copy it form the expression in the plot node for volume impedance density.
It is necessary to update your solution each time you add a variable. This is done by: Right-click your study node and select Update Solution.
These variables can now be evaluated this way:
Expand your Results node.
Right-click Derived Values, and select Global Evaluation.
In the expression field enter muscle.
Click the equal-icon. You may have to update your solution as described above.
By using this method, we get a number for the contribution from each of the model parts muscle and bone together with the contribution from the complete model. Note that it is perfectly possible to get
The expressions entered as variables may also be more complex if one wishes to expand the selectivity criteria to more than domain alone. An example is shown in the section dealing with low-quality meshes where parts of the mesh is left out. You can also use expressions such as (
The final value is the measured transfer impedance itself. This can be found in two ways: by evaluating the variable all, or by probing the voltages at the PU electrodes and divide the voltage difference by the current in the CC electrodes. This is simply the voltage difference measured on the pickup electrodes divided by the injected current. In our case, the injected current is always 1, which means that the measured voltage is the same as the transfer impedance.
To measure the voltage difference we use two boundary probes, one variable, and a Global Evaluation like this:
Expand your model node.
Right-click Definitions, and select Probes
In the Expression field, enter V (this refers to the potential in the first physics interface). The default setting is to calculate the average value of the expression or the selected boundary, which means that the probe returns the average voltage on the boundary.
Change Source Selection to Manual.
Remove all boundaries from the selection list (select all and click Remove from Selection).
Select the first PU electrode. This is PUP in figure 22, and boundary 316 in the supplied model.
Repeat 2 to 6 for the second PU electrode. This is PUN in figure 22, and boundary 315 in the supplied model.
Expand the Definitions node.
Click Variables, and enter an additional variable called transimp. Use expression bnd1-bnd2.
Right-click your study node and select Update Solution.
Expand the Results node.
Right-click Derived Values, and select Global Evaluation.
In the expression field enter transimp.
14. Click the equal-icon.
The resulting value is the same as you would get if you did the measurement using a real four-electrode measurement equipment (provided the model is accurate). An important difference is that the above method requires solving of the
As the experiment is done, it is time to have a look at the numerical results. We have found the transfer impedance contribution of muscle, bone, skin, and all materials using equation 3, and we have found the transfer impedance using voltage probes. The results are shown in table 2.
Numerical results of our experiment. The difference between
| Variable | Value | Description |
|---|---|---|
| muscle | 34+0 | transfer impedance contribution of muscle found by evaluating equation 3 for the complete model. |
| bone | 0 | transfer impedance contribution of bone found by evaluating equation 3 for the complete model. |
| skin | 0 | transfer impedance contribution of skin found by evaluating equation 3 for the complete model. |
| all | 0 | transfer impedance found by evaluating equation 3 for the complete model. |
| transimp | 1 | transfer impedance found by looking at potential at PU electrodes. |
| measured | 30 | transfer impedance measured on the real physical thigh. |
This tutorial has not looked into two-electrode impedance measurements. This is because it is very difficult to model skin impedance, effects of electrodes, and effects of polarization [1] [17].
It is possible to further reduce the model by recognizing that the skin is very thin compared to the other tissues, and not a particularly good conductor for frequencies below 1 MHz as shown in figure 7. This means that the contribution from the skin is probably negligible compared to that of the other tissues for low frequencies, and we could remove the skin and add our electrode directly to the underlying tissues. It is also possible to keep a small patch of skin under each electrode. This simplification should be done with care since the skin is becoming one of the better conductors for higher frequencies. It is advisable to verify that the effect of removing the skin is negligible before you apply this simplification on your model.
Similarly, one could argue that since bone is a poor conductor compared to the the bone marrow for low frequencies the bone marrow will not have any effect and may thus be removed from the model. In this case, the bone must be made solid by using the tool Flood Fill. This simplification is also dangerous since the relative permittivity of bone and bone marrow are starting to converge already at a few hundred Herz as shown in figure 7.
The meshing tool is adjusting the mesh coarseness according to detail level. It is possible to make the meshing tool create even coarser meshes, and thereby smaller, by adjusting the coarseness. This will reduce accuracy of the model, so care should be taken.
MPH is capable of doing a frequency sweep. This is done by replacing the frequency in the study-step nodes with an expression for a frequency range. For example, simulate from 1 kHz to 100 kHz with 4 frequencies per decade, and logarithmic frequency scale requires this expression: 10ˆrange[3,0.25,5]. Note that simulation time is 9 times that of the simulation we just have done since we specified 9 frequencies. This makes it even more important to reduce the size of the model as much as possible.
The initial values used in this tutorial had a peak value of 1500 V, while the simulation showed that values around 20,000 V would have been more appropriate since this proved to be closer to the simulated value. Furthermore, figure 26 shows that most of the voltage drop is taking place just below the CC electrodes. This means that the expression in equation 1 should have been changed to reflect this. Unfortunately, it is not always easy to know this before a simulation starts, but it may be of use for later simulations.
The tutorial has shown how to create a mesh from a MRI dataset. This mesh has been merged with electrical property data for the tissues in the mesh and used to create a model that can be simulated using the FEM software Multiphysics. A four electrode transfer impedance measurement set-up has been simulated, and a range of analyses has been done.
It is possible to create a model of a four electrode measurement system based on 3D dataset from a MRI scan or similar 3D datasets. The model allows us to explore variables that are normally unavailable to us. The data may be presented numerically and graphically.
The method may reduce the use of animal studies, and thus adhers to the three R’s of animal testing (Replace, Reduce, Refine) since some experiments can be avoided (Replace) and other experiments may be reduced and refined since we have a better understanding of what we are looking for on beforehand.