FUNDNA Software Automates Chemical Reaction Network Design for DNA Computing Circuits

A research team from the Department of Electrical and Computer Engineering at the University of Kentucky has introduced FUNDNA, an open-source design automation software tool that automatically constructs chemical reaction networks capable of executing mathematical function evaluations in biomolecular systems. Reported in the journal Frontiers in Bioengineering and Biotechnology, the platform raises the starting point of molecular circuit engineering from low-level chemical reaction specifications directly to high-level mathematical function inputs. Chemical reaction networks serve as a standard formal modeling language for molecular computation, yet engineering efficient networks that execute nontrivial mathematical operations has traditionally presented severe manual design bottlenecks. By translating user-defined single-variable functions into polynomial approximations, rearranging these expressions into modular algebraic forms, and mapping the resulting structures onto cascaded bimolecular logic gates, FUNDNA generates standardized chemical reaction networks that can be directly compiled into DNA strand-displacement architectures using existing molecular modeling tools.

Molecular computing relies on engineered biochemical substrates to execute logical and arithmetic calculations, substituting molecular concentrations for the electrical voltages typically utilized in conventional solid-state digital processors. In these biochemical systems, reaction kinetics among distinct molecular species are precisely arranged so that initial input concentrations transform deterministically into desired output concentrations over time. A substantial ecosystem of specialized software packages already exists to handle lower-level compilation and domain-level sequence modeling. Compilers such as CRNSimulator and Nuskell convert abstract chemical reaction networks into strand-displacement cascades defined by specific DNA domains and complexes, using mass-action kinetics to verify correctness and optimize kinetic topologies. Complementary suites like Visual DSD analyze domain-level strand interactions, Piperine establishes sequence designs from formal reaction networks, and NUPACK evaluates the equilibrium thermodynamics and structural dynamics of nucleic acid strands. Despite these capabilities, existing tools focus almost exclusively on processing pre-existing reaction networks rather than automating the conceptual derivation of networks that calculate user-specified functions.

FUNDNA bridges this methodology gap by operating at an architectural abstraction layer situated above the chemical reaction network level, allowing synthetic biologists to generate functional circuits without requiring extensive specialized knowledge of reaction network kinetics. The computational core of FUNDNA is grounded in an established mapping between electronic stochastic computing and molecular concentration kinetics. In electronic stochastic computing systems, continuous variables are encoded into randomized binary streams where numerical magnitude corresponds to the relative proportion of digital ones within the stream. When translated into a biochemical reaction environment, continuous values are processed via fractional encoding, in which each input, intermediate, and output variable is represented by an antagonistic pair of molecular types designated type-zero and type-one. The scalar numerical value of any variable corresponds strictly to the concentration of its type-one molecule divided by the total aggregate concentration of both type-one and type-zero species. Because this ratio represents a fractional probability bounded within the unit interval, FUNDNA assumes that all circuit inputs and outputs are normalized between zero and one.

To execute arithmetic computations within this fractional framework, FUNDNA utilizes two fundamental logic gates translated directly from stochastic computing primitives: a molecular AND gate, termed M-AND, and a molecular NAND gate, termed M-NAND. In electronic stochastic logic, a simple digital AND gate acts as an arithmetic multiplier when driven by independent randomized bit streams. Within the molecular context, an M-AND unit processes two fractional input variables represented by two distinct molecular pairs and generates an output molecular pair through a defined system of four bimolecular reactions. When type-one molecules from both inputs collide, they react to produce a type-one output molecule; all other pairwise collisions between reactant types yield type-zero output molecules. Under discrete stochastic reaction models as well as continuous mass-action kinetics, the probability of producing a type-one product is exactly equal to the product of the two input fractional concentrations. The M-NAND unit employs an identical reaction configuration among its reactant pairs but transposes the product molecular species, thereby calculating one minus the product of the two input fractions.

Before converting a target mathematical function into an arrangement of molecular gates, FUNDNA evaluates whether the candidate function complies with strict eligibility criteria dictated by the fractional encoding scheme and gate logic. Because molecular concentrations and fractional ratios cannot take negative values or exceed unity without physical normalization, target functions must satisfy specific mathematical structures. The software evaluates single-variable functions by computing a truncated Maclaurin series polynomial of user-specified degree. For algebraic rearrangement using Horner’s rule, the polynomial terms must alternate strictly in sign, and the absolute values of the coefficients must decrease monotonically as the exponent increases. These coefficient conditions guarantee that every intermediate nested term conforms to the functional form of one minus the product of two bounded fractional values. When a function fails to satisfy the alternating-sign or decreasing-magnitude requirements, Horner-based mapping cannot be applied, although the limitation reflects the structural constraints of the molecular gate library rather than an inherent algebraic restriction of Horner’s rule itself.

As an alternative rearrangement route, FUNDNA incorporates a double-NAND expansion mapping that imposes complementary constraints on the polynomial coefficients. To be eligible for double-NAND expansion, all coefficients in the Maclaurin-series polynomial must be strictly positive, and the cumulative mathematical sum of all polynomial coefficients must be less than or equal to one. These conditions ensure that the expanded algebraic expression can be recursively factored into cascaded expressions where every sub-operation corresponds to the two-input form executed directly by an M-NAND unit. When candidate functions exhibit signals that naturally fall outside the normalized zero-to-one interval, users must rescale or mathematically shift the target function prior to circuit generation, subsequently applying the inverse scaling factor to interpret simulated output concentrations. If a function fails to meet both Horner and double-NAND criteria, it cannot be processed within the current gate library without rescaling, modifying the series approximation degree, or utilizing future algorithmic expansions.

When an eligible polynomial is processed through Horner’s rule, FUNDNA systematically decomposes the nested expression into a directed graph of cascaded two-input molecular gates. In a standard Horner decomposition, the innermost parenthetical grouping matches the operational form of an M-NAND gate and is mapped directly to a single gate. Each surrounding outer parenthetical grouping contains a three-factor term combining an intermediate polynomial coefficient, the primary input variable, and the evaluated output of the preceding nested parenthesis. Because physical M-AND and M-NAND gates operate strictly on two molecular inputs, FUNDNA breaks this three-factor product into paired operations: an M-AND gate first multiplies the coefficient by the internal parenthesis output, after which a downstream M-NAND gate evaluates one minus the product of the primary input and that intermediate result. Finally, a terminating M-AND gate executes the leading constant multiplication to generate the finished fractional output. In contrast, polynomials rearranged via double-NAND expansion naturally resolve into nested pairings that map directly to successive M-NAND gates without requiring intermediate multiplier gates.

Once the gate topology is finalized, FUNDNA generates the complete chemical reaction network by substituting each constituent gate with its corresponding four-reaction chemical equations. Product molecules emerging from upstream gates serve directly as reactant inputs for downstream gates, establishing a cascaded bimolecular network where all molecular species carry explicit indices of zero or one. The complete algorithm is implemented in Python 3 and published as an open-source tool hosted in the CUT Lab GitHub repository. Through the graphical interface, users specify the single-variable function, designate the independent variable, define an evaluation point, and select the polynomial degree. Upon execution, the software renders the symbolic Maclaurin series, displays the rearranged polynomial, generates a visual graph-structure schematic of the cascaded circuit, and populates the resulting chemical reaction tables while automatically logging circuit schematics and expansions into a local assets folder.

To validate the operational capability of the automated design workflow, the authors systematically evaluated fourteen built-in function-entry templates across several mathematical families, including trigonometric expressions like sine, cosine, tangent, secant, cosecant, and cotangent; hyperbolic forms such as hyperbolic sine, cosine, tangent, secant, cosecant, and cotangent; as well as standard logarithmic and exponential templates. Beyond these built-in templates, users can enter arbitrary polynomial expressions directly into the interface provided they satisfy the coefficient requirements. From this functional suite, four representative functions meeting the mathematical eligibility constraints were selected for complete DNA-level kinetic simulations: the decaying exponential function, the hyperbolic tangent function, a quadratic polynomial defined as one-tenth of x squared plus one-fifteenth of x plus four-fifths, and a shifted trigonometric expression defined as the secant of x minus one. The quadratic polynomial and shifted secant function served as representative models for double-NAND mapping, while the exponential and hyperbolic tangent functions were processed using Horner’s rule.

The chemical reaction networks produced by FUNDNA for these four test cases ranged from sixteen reactions involving eight molecular pairs for the quadratic function up to thirty-two reactions involving thirteen molecular pairs for the decaying exponential function. To demonstrate that the resulting networks could be compiled into physical biochemical architectures, the authors used CRNSimulator to translate the bimolecular reactions into DNA strand-displacement reaction systems based on the formal mass-action kinetics framework established by Soloveichik and colleagues. In these kinetic simulations, the initial concentration of auxiliary DNA complexes was configured to ten to the negative fifth power, the maximum strand-displacement rate constant was specified at ten to the sixth per molar per second, and the baseline reaction rate constant was maintained at fifty-five thousand six hundred per molar per second. Concentrations of output molecular pairs were simulated dynamically across arbitrary time units to evaluate how rapidly the fractional ratio stabilized at its final value.

Quantitative simulation outcomes evaluated at eleven equally spaced points across the unit interval demonstrated that the DNA-level strand-displacement implementations closely tracked theoretical mathematical values. After a simulation duration of one hundred arbitrary time units, the computed fractional outputs achieved high numerical precision across all test models. For the exponential function, the computed mean square error across the unit interval was one point eight eight seven times ten to the negative seventh power. The hyperbolic tangent network yielded a mean square error of two point six one two times ten to the negative fifth power, while the quadratic polynomial achieved a mean square error of two point nine eight zero times ten to the negative twelfth power. The shifted secant function exhibited the highest deviation, reaching a mean square error of three point five one nine times ten to the negative fourth power, maintaining the overall error across all four test systems strictly below three point five two times ten to the negative fourth power.

The authors analyzed the origins of the remaining numerical discrepancies, noting that total error stems from a combination of the truncation of infinite Maclaurin series into finite polynomials, structural approximations introduced when mapping abstract reaction equations to physical strand-displacement reactions, and finite simulation durations. The largest observed deviations occurred near the upper boundary of the unit interval for functions characterized by steeper growth curves, where higher-order polynomial terms would be required for greater precision. The researchers emphasized that the current study provides computational and simulation-based validation rather than in vitro experimental execution, and that current tool operations remain limited to single-variable functions meeting strict coefficient constraints. Planned development includes integrating FUNDNA with downstream sequence-design compilers to establish an end-to-end automated synthesis pipeline, implementing automated mathematical rescaling for arbitrary signal domains, and incorporating expanded algebraic rearrangement strategies capable of accommodating broader classes of biomolecular functions.

Subject of Research: Automated chemical reaction network design for biomolecular function computing

Article Title: FUNDNA software for computing mathematical functions by DNA

Article References: Huse, J., Salehi, S. A., Ahmadi, M., & Tavasoli Naeini, M. (2026). FUNDNA software for computing mathematical functions by DNA. Frontiers in Bioengineering and Biotechnology, 14, Article 1883168. https://doi.org/10.3389/fbioe.2026.1883168

Image Credits: AI Generated

DOI: 10.3389/fbioe.2026.1883168

Keywords: FUNDNA, chemical reaction networks, molecular computing, DNA strand displacement, stochastic computing, fractional encoding, synthetic biology, Maclaurin series, Horner's rule, double-NAND expansion, CRNSimulator, biomolecular circuits, mass-action kinetics, logic gates, biochemical computation, design automation, computational biology