ENH: monopole source model

[harmening]

Enhancement: Adding the monopole source model.

• Why would you need that (since the dipole is a widely accepted approximation for cortical sources in EEG)?
  With monopoles, more complex source models can be constructed like tripoles for muscular sources or any arbitrary source configuration. The EMG of muscular sources contributes in EEG as artifacts and due to the proximity of sources and sensors is potentially better approximated by a near field tripolar solution (further details: HArtMuT).
  Thus, OpenMEEG users would not only be able to simulate cortical activity but also EEG artifacts like muscular activity more realistically.
• Key changes:
  OpenMEEG/include/monopole.h, add monopole class.
  OpenMEEG/include/analytics.h, add potential derivative.
  OpenMEEG/src/assembleSourceMat.cpp, add MonopSoureMat (analogously to DipSourceMat).
  OpenMEEG/src/operators.cpp, add operatorMonopolePot and operatorMonopolePotDer (analogously to operatorDipolePot and -Der).
  apps/assemble.cpp, add -MonopSourceMat (-MSM) (analogously to -DipSourceMat (-DSM)).
  wrapping/python/openmeeg/__init__.py, add MonopSourceMat (analogously to DipSourceMat).

[papadop]

Thank you for this contribution.

I will try to review it in the coming weeks, but it will have probably to wait a bit because I'm mostly unavailable in the next 2/3 weeks (I'll try to do my best).
But I already have one request, can you add some non-regression tests on the top of the already quite complete submission ?

Thank's again.

[papadop]

Could you give a write access to your clone of OpenMEEG, so that I can push some corrections directly on your branch. Alternatively, I can clone your branch and give you access write. What do you prefer.

[larsoner]

Or alternatively @harmening you could click the "give edit access to maintainers" button on the PR. But it's enabled by default I think so unless you un-checked it when opening the PR you should be able to push already @papadop ...

[larsoner]

@papadop did you try pushing? I think it should already work:

[papadop]

Yes I tried a simple correction and it failed, that's why I'm aslking.
Unless of course I made a mistake, I'm in a conference so my brain is shared....

[papadop]

OK. I just retried. It works....

[papadop]

OK. With my very small change, most builders pass. I will do a code review in the coming days/weeks. If, in the meantime, you can contribute the tests I was referring to, that would be nice....

Thank's again....

[harmening]

Hey, thanks for the fast reply. I added a simple test as in the dsm case, please let me know if you're satisfied or if you'd like to have sth more (sophisticated). Also haven't added any python tests for the msm case yet, could do so if you wish

[papadop]

Hello. Thank you again for your contribution.

I did a code review of your proposal and pushed some changes. I think we are mostly good, but I have a final set of questions for you....

• It seems you re-used the core of analyticDipPotDer for computing the P1term of the function. I made some changes to factorize that code between analyticDipPotDer and analyticMonopolePotDer. As I had a rough time understanding what did that code, I wonder whether you analyzed it on your side or whether you just used it blindly ?
  In particular, I discovered the hard way that the code was re-computing the barycentric coordinates of the point passed in parameter, which indeed are the values of the P1 function at that point. I ended up providing another (slightly faster but above everything clearer) of that code in the class BarycentricCoordinates.
• With respect to the previous point, I'm a little bit worried by a (very old) comment in analyticDipPotDer which questions about the sign of the returned value ("// RK: why - sign ?") and indeed I do not know why there is this sign.
  I see that you also accounted for such a sign change (because otherwise the EMpart should have had a minus sign in the code), so I wonder whether you understood why there is this sign change (probably handled at a higher level of the code, I did not investigate).
• Finally, the test you added is fine but some files must be added to allow sibling tests which are currently failing because of missing monopole files. I will handle that. This is simple and I will do that in a minute.
• On the other hand, it would be nice to compare the results of the MSM options to some groundtruth data, so we should add such comparison tests with stored results, which brings me to my last question. How have you validated this code ? Do you have analytical results on some simple model (spheres) to which we can compare the obtained results ? If not we can probably simulate a dipole with two very close monopoles and check that DSM and MSM provide similar results.... What are your thoughts on be best way to add such validation tests ?

[papadop]

Missing .monop files that prevented exiting tests to pass are now added.
Now mostly remains the question of adding comparison tests and the PR can be merged...

[papadop]

I agree with this - very old - comment.. I saw that in the analyticMonopolePotDer follows this sign change (at line 246 there should be a minus sign - per the mathematical formula - which is not there) . Did you do it to match the sign here, or have you validated by some other mean (comparison with some ground truth, or analysis of the code at the upper level) that this sign change is needed.

[DanielMiklody]

I do not really know where the sign question comes from because I do not know your derivations but I am sure it is correct that the monopole has opposite sign in this sense here. The monopole derivative is -1/4pi q@r/|r|^3 while the dipole is 1/4pi (q/|r|^3-3q@r/|r|^5*r)

[papadop]

I changed all of this code, because honestly I was not able to understand what it did from the formula. And there was no comment to explain it. The new version is explained and costs less (both in space, which is not very crucial, but also in time). I verified that it gave strictly the same results. A further improvement would be to completely eliminate this calculations. Indeed, we compute the points using barycentric coordinates and then pass the points to these functions... to re-compute the barycentric coordinates. It would be simpler to pass directly the barycentric coodinates and to recompute the point in the analytic functions. I will do this in another PR.

[papadop]

Comment for the line below. You provided the Head1/Head1.msm, which is most probably the expected result of the test. Can we add some comparison test that checks that the output of the test is indeed that one.
I'm also interested on how this file was obtained. If this is by running the OpenMEEG code, we will just check that the code always provide the same result. If the values were obtained with another numerical model, then we might check that the obtained values are indeed those that are expected. I suspect that we are in the former case... See the comments on the PR, for more suggestions on tests to add.

And obviously, we nedd that for Hed2. Head3, HeadMN*, HeadNN{a,b,c}*, ....

[papadop]

Can we also add comparison tests. See OPENMEEG_COMPARISON_TEST.

[harmening]

Thanks for you detailed code review!
I just added the expected msm results for the remaining test heads. Regarding your questions:

• We already developed the monopole potential in early 2020 and I remember we felt quite confident with understanding what we were doing back then. However, since then the architecture changed, such that for this PR I rather copied the code from the dipole potential blindly to fit the new structure, rerun the HArtMuT artifact simulations and compared the results with those of our original implementation.
• Looks good to me.
• I'm also confused with the signs and agree that mathematically there should be a minus sign. Maybe there's a minus sign somewhere else? Sth like normals orientated the other way round or so? I mean the signs in the final potential at the electrodes match the dipole orientations as one would expect.
• Yes, I accounted for this sign change in analyticMonopPotDer. Mathematically, I would have expected a plus sign (added this now as a comment).
• I indeed simulated a dipole with two very close monopoles within the HArtMuT head and checked if the summed MSMs produce similar results to the DSM. So for example the dipole [0,0,0, 0.000001*2,0,0] produced the same field as the summed monopoles [0.000001,0,0, 1] and [-0.000001,0,0, -1] up to around 12 decimals (depending on the complexity of the model).
  I'll try to translate such tests for your openmeeg test heads (also Head1, Head2, Head3, HeadMN*, HeadNN{a,b,c}* right?), but I might need some help with where to put them (I assume in tests/ComparisonTest_caller.cmake, maybe line 83?).

[agramfort]

could we have a unit test comparing the solution to some ground truth values obtained in a sphere model or with some alternative reference solver?

…

Message ID: ***@***.***>

[DanielMiklody]

could we have a unit test comparing the solution to some ground truth values obtained in a sphere model or with some alternative reference solver?
…
Message ID: @.***>

We will try to implement this. I would go for a sphere model with a single source. Do you think single or multisphere. Single is of course much easier.

[papadop]

You can add multiple tests. Sphere is fine. Both single shell ans multiple shells would be great.... If you need ground truth for multiple shells, I can probably help.

[DanielMiklody]

You can add multiple tests. Sphere is fine. Both single shell ans multiple shells would be great.... If you need ground truth for multiple shells, I can probably help.

Help wanted, I never actually simulated any spheres and don't know any solver not built for dipoles. I can work it out myself but why if you offer help? ;D

[papadop]

You can add multiple tests. Sphere is fine. Both single shell ans multiple shells would be great.... If you need ground truth for multiple shells, I can probably help.

Help wanted, I never actually simulated any spheres and don't know any solver not built for dipoles. I can work it out myself but why if you offer help? ;D

I have some software that is able to extend the solution in infinite homogenous medium to multiple spheres. I just need to investigate what is the spherical harmonics decomposition of the monpole potential....

[DanielMiklody]

I have been thinking about this today a little bit:

• One problem in general is that a single monopole as a current source is uneralistic in a sphere (or anything else) surrounded by non-conductive air. The solution however works so far. If you superimpose the solutions for two monopoles close to each other you get a dipole. I guess this is rather a theoretic problem that you externally need to guarantee that the sources always sum to 1. We have fulfilled that 'manually' in our paper. I'm just not getting yet, how the boundary conditions can be solved with unrealistic assumptions. They then guarantee that no currents leave the scalp but where do they go? 🤣
• The easiest analytical testcase would probably be a source in the center of the sphere. Like this, we have the infinite homogeneous potential at the surfaces, no need for spherical harmonic decomposition. Or am I wrong?
• In the other case, you need to care about several things: the spherical harmonics decomposition has to go down to order n=0 not 1 as in dipoles. As far as I could derive and understand it further: The homogeneous solution for the innermost sphere is actually $\frac{I}{4\pi\sigma_1 r} \sum \dots$ instead of $\frac{Id}{4\pi\sigma_1 r^2} \sum \dots$. The rest should stay the same, so all you have to do is change that prefactor of the sum in your solution formula. I have some preliminary derivations for that that I could share.

[DanielMiklody]

Hey, @papadop, are you working on a reference solution or shall I also? I got that you had something that should be easy to adopt so I was waiting but could of course also proceed implementing something. Thanks anways :)

[DanielMiklody]

we found a solution: this toolbox models mulitlayered spheres for monopoles. https://github.com/femeg/FEMEG/blob/master/femeg_anisphere.m
We have verified already that our solutions and those correlate with over 0.99.
As the implementation is in matlab, we would now do the following:

• save the leadfield from the analytic solution as a binary file
• then we need a test, that runs the same using om and compares the solution.

Are you fine with that approach? What error measures would you want? Do you have a template test that does something similar we can use to adopt? We need to run several commands to build the leadfield (om_assemble for headmodel, msm and sensor interpolation).

[harmening]

Hey @papadop, what do you think about our approach of storing the analytic monopole solution in binary files for the different heads and testing the openmeeg code against them? Or do you have any concerns? I was also thinking about implementing the analytic solution here, too, but I think that's a bit overshooting the mark, isn't it?

[harmening]

@papadop @larsoner @agramfort, sorry to bother you again, but it feels like we're stuck at this point. I would still love to see our changes merged to finally enable the community to use HArtMuTs simulation of EEG artifacts. We have pending PRs prepared for fieldtrip and eeglab's dipfit once MSM is merged.

Please let me summarize, what we have here so far:

• working code to compute the field of monopoles (similar to dipoles)
• working tests for all test heads of data folder (as done for dipoles)

We were discussing about adding further tests, like testing against analytic solutions, solutions by some alternative reference solver, or testing the solutions of dipoles to those of two very close monopoles. Since we could not proceed any further on how to implement further tests for almost one year now, what do you think about merging this branch finally and opening up a new thread/issue for the implementation of further testing. I agree that the discussed tests make totally sense and I am also willing to help implementing that (also because we already compared monopoles in spherical models against the analytic solution as provided by FEMEG), but I need help in constructing such a test within this framework here.

I am very much looking forward to your response! Would be really really great if we could somehow move forward :-)

[larsoner]

I have no objection, hopefully @papadop or @agramfort can weigh in

In the meantime I'll merge main into this PR to see if CIs are happy