cellmlmanip and chaste_codegen: automatic CellML to C++ code generation with fixes for singularities and automatically generated Jacobians

Hundreds of different mathematical models have been proposed for describing electrophysiology of various cell types. These models are quite complex (nonlinear systems of typically tens of ODEs and sometimes hundreds of parameters) and software packages such as the Cancer, Heart and Soft Tissue Environment (Chaste) C++ library have been designed to run simulations with these models in isolation or coupled to form a tissue simulation. The complexity of many of these models makes sharing and translating them to new simulation environments difficult. CellML is an XML format that offers a widely-adopted solution to this problem. This paper specifically describes the capabilities of two new Python tools: the cellmlmanip library for reading and manipulating CellML models; and chaste_codegen, a CellML to C++ converter. These tools provide a Python 3 replacement for a previous Python 2 tool (called PyCML) and they also provide additional new features that this paper describes. Most notably, they can generate analytic Jacobians without the use of proprietary software, and also find singularities occurring in equations and automatically generate and apply linear approximations to prevent numerical problems at these points.


Introduction
Within the area of electrophysiology, there are hundreds of mathematical models describing biological behaviour. Many of these models are complex systems of tens of ordinary differential equations (ODEs) with hundreds of parameters, making translation into different simulation software time-consuming and prone to transcription errors. This also makes sharing models between different tools and application areas difficult. CellML 1 addresses this problem by offering a way to describe mathematical models in an XML-based format, independent of the choice of programming language or tools used to simulate or analyse the models. It was originally created with the Physiome Project in mind, and a large repository of well over a hundred CellML electrophysiology models is available on the Physiome Model Repository 2,3 (PMR, https://models.physiomeproject.org/cellml). CellML sees continued development by the research user base 4 and several tools are available to support modelling using CellML models.
CellML models can be imported into various simulation tools such as the Cancer Heart and Soft Tissue Environment (Chaste) 5 , OpenCOR 6 , Myokit 7 , and model comparison tools such as the Cardiac Electrophysiology Web Lab 8,9 . This paper describes the development of cellmlmanip and chaste_codegen -tools that together provide a Python 3 CellML code generator for CellML model import into Chaste C++ files, replacing PyCML 10 , written in Python2. The new tools provide several useful features that we have described previously 11 (in the context of PyCML), most notably: • Automatic units conversion (electrophysiology models in Chaste always use milliVolts, milliseconds, microAmps per square centimetre for membrane currents).
• Generation of C++ code for either inbuilt Chaste ODE solvers or the Sundials CVODE library.
• Automatic generation of lookup tables for faster evaluation of voltage-dependent functions.
• Automatic generation of code for solving using Rush-Larsen style schemes.
cellmlmanip and chaste_codegen are an all-new implementation which also include: i) generation of Jacobian matrices analytically; and ii) fixes for numerical singularities. We will focus on these new features in the rest of the article.

Methods
Implementation CellML is a model definition language, as such it describes the model and its equations but does not describe how the equations should be solved or how experiments are run. The parsing of CellML files into sets of symbolic/algebraic equations (in SymPy 12 format) is handled by the cellmlmanip library. chaste_codegen takes these equations and generates C++ code using Jinja2 13 templates for compilation within Chaste. cellmlmanip is available as a Python3 library and chaste_codegen is available as a standalone command line tool for Python 3. As Python based tools, they are platform independent, and have been tested on Windows, Linux and Mac. chaste_codegen, cellmlmanip and their stack of dependencies are all free and open source. chaste_codegen (with its dependency cellmlmanip) has been integrated into Chaste v2021.1 onwards.
cellmlmanip is a flexible component that can read CellML and enable it to be used for a variety of purposes, such as translating CellML models into other formats, or into code for various simulation packages. A reason for separating this functionality into separate tools is that we want to enable parsing and checking CellML files without the need to generate Chaste code, which is useful for other projects, and indeed the next version of the Cardiac Electrophysiology Web Lab will use cellmlmanip but not chaste_codegen. Separation of the

Amendments from Version 1
The title and text have been amended to describe the roles of both the new cellmlmanip and chaste_codegen tools more clearly, and to address comments from the reviewers. The introduction contains more information on the other capabilities of the tools. There is more information on alternatives for the singularity fixing, and detail on the use of the ontology to annotate CellML models before using with chaste_codegen.
Any further responses from the reviewers can be found at the end of the article REVISED parser and simulator has the advantage of creating a resource for the CellML community, which will ultimately be much easier to maintain than writing a bespoke CellML parser for each application.
Key features of cellmlmanip are: using SymPy to represent mathematics (making its full suite of algebraic manipulation capabilities available), tracking physical units and performing conversions as needed, and managing and querying Resource Description Framework (RDF) metadata annotations on models. cellmlmanip currently supports CellML version 1.0, but it could easily be adapted to support CellML version 2.0 4 . However, because of the ongoing development of libCellML 4 there are currently no plans to do this as libCellML is a more general propose CellML 2.0 library and should offer the ability to read and manipulate CellML 2.0 models, and also write adjusted models back to CellML again, which is not something cellmlmanip aims to do. We hope to replace cellmlmanip with libCellML at some point in the future when it has all the features we need. Indeed, the singularity fixing code could migrate to libCellML in future.
SymPy is a python library for symbolic mathematics 12 that offers several convenient features. Most notably, it provides us with the ability to calculate Jacobians algebraically, to recognise patterns, rewrite equations, and to extract common terms in a set of equations. It also comes with a convenient printing mechanism, separating the mathematics from their representation.
Jinja2 is a templating language for Python, modelled after Django's templates. Using Jinja2 allows us to separate the logic from the code output, which allows generating code for a number of different solvers, as described in detail in Cooper et al. 11 . Jinja2 templates should allow easy adaptation to export code in other programming languages for other cardiac electrophysiology simulators (we are planning to use the same approach for python code generation within the Cardiac Electrophysiology Web Lab). The latest release of chaste_codegen also adds the ability to generate code in LabVIEW 14 format by adding a new Jinja2 template, highlighting the flexibility of this templating approach.
Analytic Jacobians Within Chaste one of the main solver types available is CVODE, which performs well for the stiff systems that feature in many electrophysiology models. CVODE can be sped up if the user provides a method to return the 'analytic' (algebraically derived and exact) Jacobian matrix for the ODE system, with entries defined as the partial derivative of each ODE's right-hand-side with respect to each state variable. If this method is missing CVODE derives an approximation based on finite differences 15 . To supply analytic Jacobians PyCML required them to be pre-computed by the proprietary Maple software. In chaste_codegen we have integrated calculation of a Jacobian matrix at runtime, by making use of the open source SymPy library.
We compared previously existing code, using Maple's differentiation, with the SymPy generated Jacobians. The validity of the generated Jacobians was assessed by comparing numerical results for 46 different models, all results were identical or within very small tolerances.

Singularities
Many electrophysiology models use a formulation for ion currents based on the Goldman-Hodgkin-Katz (GHK) flux equation 16 , or feature equations with similar structure. Unfortunately such equations can introduce singularities into a model. Here, we define singularities as situations where an equation tends to 0/0 close to a particular value of a model variable (usually membrane voltage, V).
In general terms, these GHK-style functions of voltage take the form: where B is any constant, f(V) is any function of V (that may also be simply a constant), and v 0 is the voltage at which we hit the singularity. Following Johnstone 17 , we can simplify this notation by defining a new function and use the substitution to leave Equation (1) with a non-dimensional fraction term that encapsulates the singularity, such that For convenience we also define We can now see more clearly that at as U → 0 there is a singularity as both the top and bottom of the fraction in Equation (5) tend to zero (e U → 1, so e U -1 tends to zero).
It is important to point out that mathematically the value of the equation remains well defined; no physics in the model is breaking down as we get close to 0/0. That is, there is a finite value to which the expression evaluates which we will tend towards as we get closer to the singularity (a limit), and this can be found analytically using approaches such as L'Hôpital's rule, as we will discuss later. But when evaluating such equations numerically, it is possible to reach voltages close to (or at) the singularity where the numerical evaluation is not just inaccurate but often tends to ±∞.
As an example, Figure 1 shows the background calcium current I Ca,b in the Davies et al. 18 model of a canine cardiomyocyte, as a function of trans-membrane voltage (V). The figure shows unstable, asymptotically-increasing oscillatory behaviour close to V = 0.
By plotting such graphs for a number of cases, Johnstone 17 found that the range in which instability occurs numerically is within approximately 10 -7 of U = 0 when using double precision (64 bits to represent a floating point number in computer memory).
Problems with these singularities are most apparent when running voltage-clamp experiments, when the voltage V is clamped exactly at a singularity voltage v 0 . Given the infinite number of choices for voltage clamps and parameters within the models, one might expect this to be unusual. However, this situation is very common.  At, and close to, V = 0mV we hit a singularity such that the computation attempts to evaluate 0/0 and answers can tend to ±∞, seen here as apparently vertical lines at V = 0 mV. Right: a zoomed-in view around the singularity. The red line represents our applied linear approximation, evaluated in the voltage interval [-1.336 × 10 -6 ≤ V ≤ 1.336 × 10 -6 ] (or equivalently -10 -7 ≤ U ≤ + 10 -7 ). The black line represents a standard fix often manually applied to models, which leads to the discontinuity between blue and black lines at U = ±10 -7 .
the (V -v 0 ) terms often feature round numbers for v 0 . As an example, equations of the form of Equation (1) with both (V + 10) and (V + 25) appeared in the original Hodgkin and Huxley 19 model, meaning that even this fundamental model requires modification to do voltage clamps to -10 or -25 mV (N.B. the voltage in the original paper is defined relative to resting potential rather than extracellular potential, and some CellML implementations have updated these numbers to the modern convention using a similar process to that described in Brown 20 ).
Apart from voltage clamp experiments and singularities at 'round numbers', any singularity within normal physiological voltage ranges will be crossed during the simulation of an action potential. Frustratingly, the more accurately you solve the ODEs, the more likely you are to hit these singularities. This is because a simple fixed-timestep solver is very unlikely to hit the narrow region of voltage which causes a problem, but using an adaptive timestep solver with automatic error correction (like CVODE 21 ) can detect large gradients or sudden changes, introduce more timesteps to refine the solution around these points, 'home in' on the singularities and crash. This is definitely not to say that the adaptive solvers are unsuitable for these models, quite the opposite -electrophysiology models frequently have stiff systems with 'fast sodium current' upstrokes and other slower processes, and adaptive ODE solvers can provide huge improvements in speed and accuracy. So rather, enabling reliable use of adaptive solvers is an added incentive to avoid the numerical problems associated with singularities.

Approaches to fix singularities
It is worth noting that the use of some computational optimisations such as using the C++ function 'expm1' (where expm1(U) = e U -1) do help by increasing precision for small values of U, but they do not alleviate the problem entirely.
There are a number of ways the singularity can be removed completely. When inspecting existing models, we found many CellML files where L'Hôpital's rule had already been applied manually when numerical problems had been encountered. In the notation we introduced above it is simple to see how L'Hôpital's rule applies So the 'fixes' in most CellML files apply this constant limit value, g(U) = 1, across a region close to the singularity (by using a 'piecewise' statement in CellML). Substituting the L'Hôpital limit into Equation (4), this fix is applied as Where ε is a small range, which varies depending on the CellML file but translates to a small region of voltage. In our notation this is equivalent to A Taylor expansion around the singularity gives 2 ( ) 1 2 12 U U g U ≈ − + + … and so as a refinement of Equation (7), Johnstone 17 suggested using the first two terms of the Taylor expansion rather than just the first term, giving a more accurate approximation close to the singularity: We have taken the approach above and implemented it in cellmlmanip with a small modification. Rather than using the above linear expression we decided to simply 'draw a line' between the values of C(V) evaluated using Equation (1) at the ε bounds of our region, and interpolate from this. This was not actually any simpler to implement, as we still identify U in order to pick bounds over which to apply the linear approximation. But a benefit of this approach over using Equation (8) is that we get at least C 0 continuity as we leave the region, even if some curvature in C(V) is apparent, whereas there can still be (typically very small) discontinuities in C(V) as we transition between the cases in Equation (8). Note that any discontinuities in Equation (8) are, by definition, much smaller than the discontinuities that appear using Equation (7), as we can see in Figure 1. This continuity can be important to avoid problems in numerical solutions of ODEs, but in practice Equation (8) and the approach we have taken are almost indistinguishable given the size of ε.
The algorithm we implemented within cellmlmanip has two parts. To identify whether an equation has a singularity and apply a fix, the algorithm is: 1. Recursively check within each equation for singularities, term by term.
2. Skip (sub-)equations that are piecewise statements, we assume that if these have a singularity it has a manual fix applied.
4. Solve for U = 0 to find the singularity point in terms of model variables (usually Voltage, V).
5. Introduce a piecewise statement to replace the original expression within -10 -7 ≤ U ≤ +10 -7 (the fact this is now a nondimensional range means that the same fixed range appears to work well for all singularities we have found; it is translated back into voltage ranges of different widths at code generation time).
6. Within this range we use linear interpolation between values evaluated using the original expression at the boundaries (U = ±10 −7 ). 3. We start at the bottom level of the graph, look for singularities and if none are found using the procedure above we progress to the next level of the graph. If a singularity is found we introduce a fix at that node of the graph.

Testing
To test the process we searched a set of CellML models (all those annotated for use with Chaste and the WebLab and available at https://github.com/chaste/cellml 22 ) for any piecewise elements. We then manually identified the subset of these being used to fix singularities. We removed these fixes from the CellML files then verified that our code would indeed find and fix these singularities automatically. The result of this exercise is shown in Table 1. The table shows the number of previously hard-coded fixes and the number of  extra fixes detected for the range of CellML files. In short, all 106 singularities with previously hard-coded fixes were identified, and the code found 263 new singularities that had never been 'fixed' within the CellML files. An overview of all 369 fixes, each with a similar plot to Figure 1, is available in the 'assets' branch of the chaste_codegen GitHub page 23 .
Note that the figure of 106 singularities includes approximately 20 hard-coded fixes we applied to our subset of CellML models using Equation (8) after hitting singularities in previous work, and these do not feature in the main CellML repository (PMR) which contains 80 to 90 hard-coded fixes. There are a number of advantages in removing hard-coded singularity fixes from the CellML files: we can make more accurate linear rather than constant fixes (as shown in Figure 1); others could choose how to treat the singularities and adapt our code to use other methods if they wish; and finally the CellML model is simplified to represent the equations and not how they should be solved. A disadvantage in removing hard-coded singularity fixes from the main CellML repository (PMR) is that all translation code would need to adopt an approach similar to the one we have outlined to avoid hitting singularities, whereas currently around 80 to 90 hard-coded fixes will appear in any generated code without code generation tools needing to do any special treatment. Dealing with singularities at code generation time is more future-proof for when people add new models, rather than hard-coding the fixes we applied here into the CellML files in the Physiome Model Repository (PMR) and then having to periodically repeat this exercise for any new models that have been added. On balance, we think that this automated approach at code generation time is the preferable route and hope it has been described so that others can reproduce it in their own code generation software, or re-use components of our open source implementation.
In the remainder of this article we discuss the practicalities of using chaste_codegen. The main limitation of this approach is that it is only capable of fixing singularities caused by GHK flux-like equations. We have decided to limit to these as they are most common type in the cardiac modelling area, and appear to cover all singularities we could find in the CellML electrophysiology models in Table 1; whilst SymPy can detect singularities more generally, this is much more computationally expensive than our 'pattern matching' approach, and so covering all possible singularities would be computationally intensive, and does not appear to be necessary for electrophysiology models.

Operation
There are two main ways to use chaste_codegen. It can be used as part of the Chaste 5 build process, from the 2021.1 release onward. It can also be used as a standalone command line tool or Python library. The minimal system requirements for chaste_codegen as a standalone command-line tool are: • python3 (3.5+), tested on Windows 10, Ubuntu Linux 18.04 and MacOS.
• python3 venv (or other python virtual environment) is recommended, to ensure the right versions of dependencies are available. However chaste_codegen will still work without a virtual environment. Python3_venv is required for use within Chaste and is usually bundled with a python installation.
For use integrated into the Chaste 5 build process, follow the regular guides 24 on installing and building Chaste. As part of the Chaste installation process chaste_codegen will be installed in a virtual environment and all CellML files in the source will be converted using the appropriate settings.
To install as a stand-alone command line tool, run pip install chaste_codegen. You may want to create a python virtual environment (venv) first. The basic usage is: chaste_codegen cellml_file where cellml_file is the CellML file to be converted. To get a detailed overview of the various options run the command chaste_codegen -h. Before using CellML files with chaste_codegen, they are required to be annotated with metadata in RDF format, following the Web Lab Ontology (https://github.com/ModellingWe-bLab/ontologies). The annotations enable variables such as stimulus current, time and voltage to be identified and then converted to consistent units. See the Chaste wiki 25 for more details.

Use cases
In this section we will briefly show a number of common conversions in action. For a more detailed guide see the 'Code generation from CellML' section in the Chaste guides 24 . Due to the size of both the import and generated code, however, we will mention only key snippets and refer to the full files available in the assets branch of the chaste_codegen GitHub page 23 .
In the following examples a CellML 1.0 file for the Hodgkin-Huxley model will be converted into code for a number of different ODE solvers. The CellML file consists of a number of components (membrane, sodium channel, sodium channel m gate, sodium channel h gate, potassium channel, potassium channel n gate and leakage current) and within these components variables are defined. The model also includes and links between components and unit definitions.
CellML files may include metadata through the use of RDF, the Resource Description Framework. chaste_codegen makes use of these annotations when generating C++ source code for Chaste, some of which are optional. In particular, chaste_codegen needs to know which variables represent voltage and stimulus_current, in order to link the models into the mono/bi-domain equations.
Below an example of a variable called V tagged as voltage is shown.
• GetIIonic calculates total ionic current at the present time (for use in tissue simulations).
• EvaluateYDerivatives calculates the derivatives of the state values when provided with their current values, defining the ODEs of the model.
• ComputeDerivedQuantities gives a way to calculate the value of quantities that are derived directly from state variables, for example currents such as "the fast sodium current".
• OdeSystemInformation::Initialise gives a way to retrieve information about the model such as name, free variable (usually time), state variable, modifiable parameters and named derived quantities.

CVODE
The following command generates code for the CVODE solver which has its own vector class.

chaste_codegen -cvode -use-analytic-jacobian ModelName.cellml
This generates .cpp and .hpp files with the same name as before. The generated class now inherits from AbstractCvodeCell, containing the same methods but with SUNDIALS' vector class for use directly with CVODE. It also has a method EvaluateAnalyticJacobian in which the analytic Jacobian is defined, to be used by CVODE.

Backward Euler
The following command generates Backward Euler code.

chaste_codegen -backward-euler ModelName.cellml
This generates .cpp and .hpp files with the same name as before. The generated class inherits from AbstractBackwardEulerCardiacCell. The class does not have EvaluateYDerivatives, but instead it has: • UpdateTransmembranePotential where the voltage is updated based on the ODE for voltage.

Rush-Larsen
The following command generates code in which some state variables are updated using analytic solutions using the Rush-Larsen scheme 26 . • ComputeOneStepExceptVoltage where the linear state-variables are updated using their analytic solutions.

Summary
This paper has introduced the cellmlmanip library for reading and manipulating CellML models and the chaste_codegen software tool, designed to translate electrophysiology models from the CellML XML format into C++ code for use by the Chaste simulation package. We have shown how the tool can be used, highlighted the main different types of code it can generate for different solvers and shown a number of advanced features that chaste_codegen implements over a previous tool called PyCML. The most notable are the ability to generate analytic Jacobians and to evaluate and fix singularities in equations. We have shown that the singularity analysis works as expected with an analysis of a large number of popular models, by identifying all previously-identified singularities and finding over three times as many in total.

Contributing to development
chaste_codegen and the cellmlmanip libraries are open-source and publicly available and we welcome contributions: from questions about the current functionality to suggestions for improvements and source code contributions. It should also be relatively simple to extend the use of templates to generate code for other simulation packages in C++, Python, or other languages. Contributions are made in the first instance using GitHub issues. In order to contribute, users create a new issue in either GitHub repository or comment on an existing issue. Contributions in the form of source code can be made by issuing a pull request on either repository (ideally a new GitHub issue should be created, which can then be linked to the pull request). The pull request will trigger an automated test suite and a number of other checks for things such as code formatting and test coverage.

BSD 3-Clause License
Author contributions MH led the development of chaste_codegen, contributed to cellmlmanip development and co-wrote this paper. MC contributed to cellmlmanip development and provided code reviews and suggestions for chaste_codegen development. AUT and SMK contributed to cellmlmanip development. RHJ identified the singularities issue and provided analysis around how this could be solved. JC led the development of cellmlmanip and provided code reviews and suggestions to chaste_codegen. GRM worked on the singularity fixing algorithm, supervised the work, and co-wrote this paper.

Open Peer Review I confirm that I have read this submission and believe that I have an appropriate level of
as suggested, but couldn't export the cpp file from the CellML model. It will be useful to provide some sort of "test_codegen" command to run sample tests and report errors to help to troubleshoot.
A concern on the tool is that it is based on cellmlmanip which supports only CellML 1.0 and the author described there are no plans to support CellML 2.0. Would this mean the chaste_codegen will be outdated when CellML 2.0 is adapted widely? It would be useful to discuss this in the article, as a part of the limitation if so. Also, would chaste_codegen throw a clear warning if CellML 2.0 is used?

5.
I would suggest authors, to thoroughly check all the equations (1 -8) and ensure all components are properly defined. For example, V0 is defined as the voltage at which we hit the singularity, but there is no definition for B in equation 1. Although, it might be trivial to those working in electrophysiology modeling and references are provided, it will be useful to describe functions such as C(V), A(V) in the text to help general readers.

6.
Validity of the singularity detection by chaste_codegen was shown through the comparison of the number of singularities detected and the hardcoded fixes in previous CellML models (Table 1). Similarly, the authors described that the validity of the generated Jacobians was assessed by comparing numerical results of runs of a number of different models without providing further details. If not in detail, at least a few sentences on how many models were used and so on would be beneficial.

7.
The authors have discussed the advantages of the tool and the newly developed features, it would be useful to also discuss limitations if any. 8.
It will be also useful if the authors discuss any comparable existing approaches and implementation for automatic singularities detection and fixing. 9.

Is the rationale for developing the new software tool clearly explained? Yes
Is the description of the software tool technically sound? Yes

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes
Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Yes Reviewer Expertise: I am a mathematical modeler and one of the editorial board members of the community standard, SBML. I am familiar with XML-based formats to describe models and the libraries developed to handle them. I also lead the BioModels, a repository of mathematical models of biological systems.
I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard. We have expanded the introduction to point out the principal features of PyCML which indeed chaste_codegen also provides.
Also, it would be useful to discuss the performance of chaste_codegen when compared against PyCML. For example, it will be useful to briefly mention how well does the computational time compares in general. Were the functionalities from PyCML comparable to those in chaste_codegen excluding the new features (automated Jacobians and singularity fixing)?
We haven't done a thorough comparison of performance other than that it was fast enough and did not noticably impact the CellML conversion time during Chaste compilation (in either direction!). As the reviewer has pointed out, chaste_codegen has a number of features PyCML misses. The primary concern for users is almost certainly run time of the generated code (solving ODEs) rather than this code generation step.
The article describes in detail the rationale and solution to fix singularities in the equations (where the value of variable such as Voltage tends to become 0/0). The chaste_codegen has been already integrated into chaste. Authors could discuss how this tool or approach could be used beyond chaste and provide a few suggestions.
We have made some suggestions about using Jinja2 templates and pointed to their documentation.
I tried to run and test the standalone chaste_codegen python tool. Although I tried only for an hour or so, I couldn't get it working. I managed to get it installed in a virtual environment as suggested, but couldn't export the cpp file from the CellML model. It will be useful to provide some sort of "test_codegen" command to run sample tests and report errors to help to troubleshoot.
Running the command chaste_codegen -h gives detailed help. Having said that, chaste_codegen should convert it. Chaste_codegen has been tested against all the models in https://github.com/Chaste/cellml It would be interesting to know what errors the reviewer faced, and we would encourage them to leave an issue on the github repository. We do require the models to have metadata annotation which is not present in most models in the CellML repository. We have highlighted this in the paper.
A concern on the tool is that it is based on cellmlmanip which supports only CellML 1.0  These were generic functions or constants that could take any form to fit particular models, we have attempted to clarify this in the text.
Validity of the singularity detection by chaste_codegen was shown through the comparison of the number of singularities detected and the hardcoded fixes in previous CellML models (Table 1). Similarly, the authors described that the validity of the generated Jacobians was assessed by comparing numerical results of runs of a number of different models without providing further details. If not in detail, at least a few sentences on how many models were used and so on would be beneficial.
We added the number of tests and explained the results in more detail.
The authors have discussed the advantages of the tool and the newly developed features, it would be useful to also discuss limitations if any.
We have added a few sentences on limitations, mainly in terms of the singularity detection only working for the GHK-flux style equations we discuss here.
It will be also useful if the authors discuss any comparable existing approaches and implementation for automatic singularities detection and fixing.
We have not seen any comparable implementations for automatic singularity fixing, all existing approaches have relied on spotting these and fixing manually as far as we are aware.

David Nickerson
Auckland Bioengineering Institute, University of Auckland, Auckland, New Zealand This article very nicely lays out the rationale for the chaste_codegen tool, particularly in the evolution from the previous PyCML software used in Chaste. While clearly the emphasis is on producing code suitable for use in Chaste, the authors do suggest how the tool could be used independently. Some introduction or pointers as to how readers might be able to edit/alter the Jinja2 templates to produce non-Chaste code would be helpful.
While the performance and correctness of the analytic Jacobians generated by chaste_codegen using SymPy is very briefly mentioned, seeing some of those checks would I think be interesting to readers.
A large portion of this manuscript is describing the use of cellmlmanip to detect and fix singularities. While this is a valuable contribution the authors preface this by suggesting the replacement of cellmlmanip with the libCellML when looking at handling more recent versions of the CellML specification. Perhaps a bit more clarity regarding the future use of cellmlmanip and chaste_codegen regarding the singularity detection and fixing would be helpful.
I needed the assistance of google to learn what expm1 is, so perhaps other readers might appreciate a pointer.
Step 2 in the algorithm to identify singularities in equations seems a bit too general. While historically most cardiac electrophysiological models have been encoded in CellML in a similar manner for which this assumption holds true, it is not likely to always be true. Particularly if readers were looking to build on chaste_codegen to work with models from other domains.
The discussion regarding the removal of hard-coded singularity fixes from the original models available in the main CellML repository is great to see. And something that the authors should be encouraged to discuss with the broader CellML community.
(Minor comment: `chaste_codegen -h` suggests that the cellml_file could be a URI, but pointing to, for The resulting Jacobians from Maple and from sympy look quite different. For a few simple toy examples we manually verified them, but beyond that, the verification is based on the converted models using the sympy generated Jacobians passing the same set of numeric tests against reference results as we previously ran for the Maple-generated Jacobians. Some further information about numbers of models tested has been added.
A large portion of this manuscript is describing the use of cellmlmanip to detect and fix singularities. While this is a valuable contribution the authors preface this by suggesting the replacement of cellmlmanip with the libCellML when looking at handling more recent versions of the CellML specification. Perhaps a bit more clarity regarding the future use of cellmlmanip and chaste_codegen regarding the singularity detection and fixing would be helpful.
While this functionality is indeed in located cellmlmanip at the moment, it would be quite straightforward to move to chaste_codegen, or elsewhere such as a future version of libCellML. It does not rely in cellmlmanip specifically, but just requires an algebraic representation of the equations such as SymPy uses. We added text to explain this in the Implementation section.

I needed the assistance of google to learn what expm1 is, so perhaps other readers might appreciate a pointer.
We add an explanation in brackets when it is first encountered.
Step 2 in the algorithm to identify singularities in equations seems a bit too general. While historically most cardiac electrophysiological models have been encoded in CellML in a similar manner for which this assumption holds true, it is not likely to always be true. Particularly if readers were looking to build on chaste_codegen to work with models from other domains.
It is indeed a limitation that we will only spot these GHK-flux style equations. One nice feature is that we still spot them even when intermediate equations are involved with our recursive approach. It would be possible to do further 'pattern matching' approaches on a case-by-case basis, or to use SymPy's inbuilt methods ( https://docs.sympy.org/latest/modules/calculus/index.html#sympy.calculus.singularities.singularities ) to find singularities in almost any expression, we tried this initially but found it was a lot slower than our approach (minutes rather than less than a second for a moderately large action potential model, and we even found some examples where the method didn't terminate at all).
The discussion regarding the removal of hard-coded singularity fixes from the original models available in the main CellML repository is great to see. And something that the authors should be encouraged to discuss with the broader CellML community.
Thanks for that suggestion, we will be attending CellML workshops and COMBINE etc. to discuss these issues.