Implementation of Magnetic Dipole representation in ModelTMSCoil

[johaenns]

@jessdtate

The implementation of the coil representation through magnetic dipoles in ModelTMSCoil is very confusing (and to my understanding partially wrong).

The goal at this point is to represent a circular coil by its magnetic dipole moment to be used as input to SolveBiotSavart.

A pointcould representation of the coil is created (similar to the "Wire" representations) and each of these points is assigned with a magnetic dipole moment. The moments for each layer of the coil are assigned according to the area in between two neighboring windings.

There are two problems with this:

1. The overall strength of the magnetic dipole moment for one loop is calculated to be the area in between two loops times the current strength. To my understanding this should be the overall area that is enclosed by each loop times the current strength, the current calculation would only be correct, if the current in the two loops had opposing signs.

2. The overall magnetic dipole moment of a coil is afterwards split up into "subdipoles" on the radius of the coil, which are then used in SolveBiotSavart. I am not sure if this procedure is accurate, it seems that, the representation of the actual field generated by the coil is still pretty inaccurate.

Point 1 is relatively easy to fix, the calculation of "ringArea" in Class DipolesCoilgen in ModelGenericColiAlgorithm.cc would have to be adjusted.

Overall, the question is whether this is the behavior the user expects. The representation of a coil through elemental magnetic dipoles is pretty tedious and does not align with the other two options that this module provides. It might be easier to understand, if instead a vector representation of the electric field is generated and SolveBiotSavart is adapted accordingly.

[dcwhite]

@pip010 was the original author of the modules mentioned.

[pip010]

Hi all,

indeed the the dipole presentation of a coil is problematic, however not in the way @jessdtate summarized it here. Before committing the dipoles model, I explained to @moritzdannhauer that the module in its current form is missing the weighting factors that will supposedly deliver correct magnitude of the magnetic field, thus my reluctance to commit it. The algorithm is directly from the paper of A.Thielscher "Linking Physics with Physiology in TMS: A Sphere Field Model to Determine the Cortical Stimulation Site in TMS" , however (even in the appendix) it is missing a formula of how those weighting factors per dipole are to be derived in the general case and I failed to see any patterns and come with an math-inductive definition.

Still, the shape of the fields are correct when compared to the thin-wire solver, which BTW is validated for both the magnetic flux B and magnetic vector potential A. However, both are quite under estimated in magnitude in the case of ~300 dipoles (which according to the before-mention paper should be enough).

So my suggestion was to either ask for help from the authors of the paper or simply drop it. The third, volumetric coil model relying on current density is also not validated in any way, so I also suggest to drop it before we can demonstrate that it can deliver good approximation in magnitude from a CAD derived 3D model of a realistic coil.

Side note: I do plan to help further with the modules and might need to redesign them. For example, what is really confusing for a user IMHO (also for @moritzdannhauer ) is that the current for the B and A field are not the same. If a user wants the A field, which is to be used for FEM to derive the scalar potential, one needs to scale the current in accordance to its temporal derivative, e.g. if for B a current of 5000A goes through the coil then for A it needs to be scaled by a factor of 1sec/0.0002sec since a typical TMS pulse has a ramp-up time of ~ 0.2 milliseconds, giving us 25x10E^6 A in the quasistatic regime. Ideally we should have LRC circuits simulator with all parameters given from a particular TMS manufacturer to be filled in the GUI of the module in question.

I am open to suggestions and even voip talk to clarify further some of the points. Let me know!

[johaenns]

Hi Petar,

thanks for your answer! Besides the implementation details - when I directly compare the B-fields computed with the Magnetic Dipoles and the Thin Wire models, the shapes are kind of similar. However, this is also the case if I just place a single magnetic dipole in the coil center. I have problems following the reasoning for placing the evaluation points on the winding to perform an integration over an area - to my understanding these should be placed inside the area (as it is also the case in the old scenarios that are still in the brain stimulator tutorial I think). In the current representation, also the computational advantage of the dipole representation vanishes?

To my understanding, the basic misconception in the adaption of the Thielscher paper is the integration approach. As it is visible in their Fig. A2, the integration has to be performed over the whole area enclosed by a single winding, which also fits to the definition of the magnetic dipole moment.

Therefore, the two mentioned errors in the implementation:

1. The overall area is the whole area that is enclosed by a single winding, not the area in between two windings. This leads to an error in the magnitude computation.
2. The points at which the magnetic dipoles are placed should be inside the area enclosed by the winding for an optimal approximation, not on the winding boundary. This is comparable to a numerical integration over an area.

Besides this, I agree that the implementation is problematic since Thielscher are referencing another paper for the construction of the evaluation points and weighting factors, which also doesn't directly contain the necessary information.

As long as we are not able to compute the correct weights and evaluation points, a single dipole representation of the coil might even be the better option?

Regarding your side note - the module at the moment states that you should enter dI/dt. If I am right, for getting the B-field instead of dB/dt or the A-field instead of dA/dt I would have to enter I instead of dI/dt at this position. Is that correct?

[pip010]

Hi johaenns,

directly on your points: 1) From what I understand from your explanation, I should simply not divide by the number of elements (ring segments) ? I will try it out in my own spin of scirun4 first let you know!

this directly from the paper: "To summarize, the length of a dipole was determined by the area of the subsurface and n times the coil current, whereby n represents the number of loops surrounding the dipole." what is the subsurface in the case, the whole ring/winding ? Because: "In a first step, each wing of the coil was modeled as a circular disk and divided into 16 rings. Each ring was further divided into elements, as shown in Fig. A2. " I assumed the are under one small ring segment is the element.

2) Indeed I checked and you are correct that I might not take the correct place (xyz) of the dipole, however the main problem is in the dipole moment!? I will try submit a bug-fix pull request, together with a small mistake in the position of the layers too!

About your single dipole suggestion: I guess the main issue is the approximation of the spread of the field coming from the surface are of the coil (figure-of-8), that is, it's spatial distribution. However, this is trivial to check too.

About the GUI for B and A : for the B field you need to enter the max current Imax, where for the A you enter the dI/dt in accordance whether you have bi-phasic or mono-phasic TMS stimulator.

p.s. is this where the latest of the brain stimulator toolkit is : https://www.sci.utah.edu/cibc-software/scirun/brainstimulator.html ??? I can look into it !

[johaenns]

Hi Petar,

1. The subsurface is the area of the coil that is represented by a single dipole, the construction of the subsurfaces is shown in Fig. A2. Each subsurface is represented by a single dipole with strength (area of subsurface)current(number of windings). I think I understood the subsurface construction now, I'll give it a try to implement it.

2. The single dipole would not provide a good representation of the field close to the coil, that's correct.

[johaenns]

Ok, so after a bit more reading, this seems to be what Thielscher et al. do:

The approach of Thielscher et al. is based on Ravazzani et al., 1996, who base their work on Roth and Sato, 1992 (https://www.researchgate.net/publication/277711170_Accurate_and_efficient_formulas_for_averaging_the_magnetic_field_over_a_circular_coil). The original idea is to use the Stroud integration points and weights (Stroud, 1971, Approximate calculation of multiple integrals), in the both paper above up to 7th order = 12 points, to compute the integral over the coil surface area.

Thielscher et al. expanded these formula:

They subdivide the inner coil area and weight it by the number of windings. Therefore, they subdivide the area into rings and divide each ring into multiple segments. Further they subdivide the space in between the windings and weight it by the number of windings that are still enclosing this area. The used radii and number of segments for each ring are found in Table A1.

Unfortunately, it does not get clear from the Thielscher et al. paper, how they derive the number of rings, how many segments to use, etc. to be able to adapt this for general coils.

Since I do not have time to look deeper into this right now, I would for now implement the single dipole approach, which should at least be more correct than the current implementation.

[pip010]

Unfortunately, it does not get clear from the Thielscher et al. paper, how they derive the number of rings, how many segments to use, etc. to be able to adapt this for general coils.

Exactly the problem I faced when coming with the initial implementation! Thanks for digging, will have a look into Roth's paper.

Try with a single dipole, though I expect both the magnitude and shape to be problematic. You can take the thin-wire model as a golden standard for now.

[github-actions[bot]]

Stale issue message

[dcwhite]

@jessdtate Beta or post-beta?

[jessdtate]

I think @johaenns fixed this before he left, actually. let me double check and see if we can close this

[github-actions[bot]]

This issue is stale because it has been open 120 days with no activity. Remove stale label or comment or this will be closed in 60 days.

[dcwhite]

@jessdtate Did you check whether something was done on this one?

[jessdtate]

1805 Says it might not be fixed

[dcwhite]

OK

[jessdtate]

The single dipole for a double loop puts 2 vectors in the center of both coils, rather than the center for each of the respective coil.

[github-actions[bot]]

This issue is stale because it has been open 240 days with no activity. Remove the stale label or comment, or this will be closed in 60 days.

[github-actions[bot]]

This issue is stale because it has been open 240 days with no activity. Remove the stale label or comment, or this will be closed in 60 days.