I was using ModelingToolkit, and had a system set up, but for some reason, even in a simple example, the results from the hand-spun system and the MTK one were pretty different, and the solving was also pretty slow (understandably, since the constraints caused the final equation to be a DAE rather than an ODE).
After my exams end, my plan is to look into constructing the system on an expression level, by appending coupling terms to the equations. Since the metadata of the system (jacobian, Wfact, etc) seem not to be constructed, I think it should be possible to construct the ODE systems, and then append to them.
I’m using the FitzHugh-Nagumo model, by the way, which is simpler than most of the more commonly used ones but still reasonably accurate.