We recently designed an innovative scaffold-bioreactor unit for the bioengineering of a three-dimensional (3D) bioartificial human thyroid gland or its miniaturized replica as a part of a microfluidic chip test system. This device is based on the evidence that the 3D geometry of the intraglandular stromal/vascular scaffold (SVS; i.e., the fibrous and vascular matrix) of mammalian viscera plays a key role in guiding growth and differentiation of in vitro seeded cells. Therefore, we initiated a research program focused on computer-aided reconstruction of the 2nd to 4th order intralobar arterial network (IAN) of the human thyroid gland as a reliable surrogate for its 3D SVS, to be used as an input for rapid prototyping of a biomaterial replica. To this end, we developed a computational template that works within the Mathematica environment, giving rise to a quasi-fractal growth of the IAN distribution, constrained within an approximation of the thyroid lobe shape as a closed surface. Starting from edge detection of planar images of real human thyroid lobes acquired by in vivo real-time ultrasonography, we performed data approximation of the lobar profiles based on splines and Bezier curves, providing 3D lobar shapes as geometric boundaries for vessel growth by a diffusion-limited aggregation model. Our numerical procedures allowed for a robust connection between development of lobar arterial trees and thyroid lobe shape, led to a vascular self-similarity consistent with that of a cadaveric lobar arterial cast, and reproduced arterial vessels in a proportion not statistically different from that described for the real human thyroid gland. We conclude that our algorithmic template offers a reliable reproduction of the extremely complex IAN of the adult human thyroid lobe, potentially useful as a computational guidance for bioprinting of thyroid lobe matrix replicas. In addition, due to the simplicity and limited number of morphometrical parameters required by our system, we predict its application to the design of a number of patient-tailored human bioartificial organs and organs-on-chip, including parenchymal viscera and bones. Impact statementThe study introduces the computer simulation of the three-dimensional (3D) intrinsic vascular matrix of the human thyroid gland, offering a general concept applicable to a number of other human viscera. Indeed, it provides a flexible software tool for reproduction of a 3D surrogate of the organ's 3D stromal matrix, suitable for eventual 3D bioprinting with biomaterials, and recellularization with organ-specific stem cells/progenitors. The final expectation is the design of patient-tailored 3D organ's matrices upon clinical request.

The long-accepted thermoelectric model for electro-discharge machining is being brought into question. Several experimental facts prompt the proposal of a new theory based on the effect of gap pollution on the ignition of discharges. The first experimental proof comes from the recently reported observation of debris chains and clusters. In this view, each step of the process depends on the previous ones through a deterministic relation, even if the overall evolution is unpredictable. The paper establishes mathematical grounds for the abovementioned intuitions by setting up and solving a recursive equation for the machining energy employed at each discharge event. By means of numerical and algebraic tools, the above equation is studied and shows a chaotic evolution similar to that of the logistic map. Results reconcile the apparent paradox between deterministic nature and stochastic localization of the discharges and introduce a description of the chaotic dynamics of electro-discharge processes.

It is currently known that a number of human vascular systems have a fractal geometry. Since we have recently developed a technique to prototype single arterial branches of human soft tissue organs by additive layer manufacturing (AM), we have explored the possibility that auto-similarity in vessel branching represents a key variable for accurate computational modeling of the organ three-dimensional (3D) macro/microscopic anatomy, and its reproduction by inverse engineering. To this purpose, ramification features of the intralobar arteries of the human thyroidwere studied using injection-corrosion casts of the cadaveric gland. Vessel diameters, ramification angles, and branch lengths were measured by light microscopic, computer-aided optical metrology. Distribution of morphological variables was considered on a cumulative basis, and special focus was given to the branching laws. To reduce the bias of vascular distortion due to the pressure of intravascular resin injection, measures were made dimensionless through the use of a scaling parameter set on the vascular caliber of major afferent arteries. In addition, using high resolution microtomography (mCT Skyscan 1172, Bruker microCT) equipped with CTAn software and the Otsu algorithm for segmentation, spaces occupied by vascular branches (referred to as Volume of Interests, VOI) were selected, and their planar fractal dimension calculated. Finally, a computational simulation of the vascular tree was achieved using a mixed, stochastic/deterministic algorithm, based on diffusion limited aggregation (DLA), constrained by mean values of vascular variables. Ratios among decreasing cast calibers, ramification angles, and branch lengths, respectively, were found strictly interrelated, mCT-VOI depicted fractal dimensions, and DLA simulation led to a fractal-like organization consistent with real data morphometrics. In summary, thyroid arterial geometry reliably exhibited a degree of auto-similarity, suggesting that fractality is a key feature for computational modeling and eventual AM of 3D vascular networks of the human thyroid.

This article outlines design and implementation details of the framework for one step methods for solving ordinary dieren tial equations in Mathematica. The solver breaks up the solution into three main phases for equation processing and classication, numerical solution and processing of results. One of the distinguishing features of the framework is the hierarchical nature of the method invocation which allows for simple construction of composed integration schemes. A plug-in facility for user dened schemes is also provided. Highly accurate reference solutions can also be obtained by making use of arbitrary precision software arithmetic. Issues relating to appropriate formulation and ecien t implementation will also be discussed, together with strategies for automatic method, order and parameter selection. c 2008 European Society of Computational Methods in Sciences and Engineering

Arithmetic systems such as those based on IEEE standards currently make no attempt to track the propagation of errors. A formal error analysis, however, can be complicated and is often confined to the realm of experts in numerical analysis. In recent years, there has been a resurgence of interest in automated methods for accurately monitoring the error propagation. In this article, a floating-point system based on significance arithmetic will be described. Details of the implementation in Mathematica will be given along with examples that illustrate the design goals and differences over conventional fixed-precision floating-point systems.

Composition and splitting are useful techniques for constructing special purpose integration methods for numerically solving many types of differential equations. In this article we will review these methods and summarise the essential ingredients of an implementation that has recently been added to a framework for solving differential equations in Mathematica.

This work focuses on the derivation of composition methods for the numerical integration of ordinary differential equations, which give rise to very challenging optimization problems. Composition is a useful technique for constructing high order approximations, while conserving certain geometric properties. We survey existing composition methods and describe results of an intensive numerical search for new methods. Details of the search procedure are given along with numerical examples, which indicate that the new methods perform better than previously known methods. Some insight into the location of global minima for these problems is obtained as a result.

Explicit Runge-Kutta schemes are the methods of choice for solving nonstiff systems of ordinary differential equations at low to medium tolerances. The construction of optimal formulae has been the subject of much research. In this article, it will be shown how to construct some low order formula pairs using tools from computer algebra. Our focus will be on methods that are equipped with local error detection (for adaptivity in the step size) and with the ability to detect stiffness. It will be demonstrated how criteria governing ‘optimal’ tuning of free parameters and matching of the embedded method can be accomplished by forming a constrained optimization problem. In contrast to standard numerical optimization processes our approach finds an exact (infinite precision) global minimum. Quantitative measures will be given comparing our new methods with some established formula pairs.

We give an overview of the construction of algebraic conditions for determining the order of Runge-Kutta methods and describe a novel extension for numerically solving systems of differential equations. The new schemes, called Elementary Differential Runge-Kutta methods, include as a subset Runge-Kutta methods, Taylor series methods, Multiderivative Runge-Kutta methods. We outline how order conditions have been constructed for the new schemes using B-series and their composition and give details relating to a Mathematica implementation.

Strategies for reducing the effect of cumulative rounding errors in geometric numerical integration are outlined. The focus is, in particular, on the solution of separable Hamiltonian systems using explicit symplectic integration methods and on solving orthogonal matrix differential systems using projection. Examples are given that demonstrate the advantages of an increment formulation over the standard implementation of conventional integrators. We describe how the aforementioned special purpose integration methods have been set up in a uniform, modular and extensible framework being developed in the problem solving environment Mathematica .

This paper focuses on the solution of separable Hamiltonian systems using explicit symplectic integration methods. Strategies for reducing the effect of cumulative rounding errors are outlined and advantages over a standard formulation are demonstrated. Procedures for automatically choosing appropriate methods are also described.

A component of a new environment for the numerical solution of ordinary differential equations in Mathematica is outlined. We briefly describe how special purpose integration methods can be constructed to solve structured dynamical systems. In particular we focus on the solution of orthogonal matrix differential systems using projection. Examples are given to illustrate the advantages of a projection scheme over conventional integration methods.

Many systems of ordinary differential equations are quadratic: the derivative can be expressed as a quadratic function of the dependent variable. We demonstrate that this feature can be exploited in the numerical solution by Runge-Kutta methods, since the quadratic structure serves to decrease the number of order conditions. We discuss issues related to construction design and implementation and present a number of new methods of Runge-Kutta and Runge-Kutta-Nyström type that display superior behaviour when applied to quadratic ordinary differential equations.

Much recent work has indicated that considerable benefit arises from the use of symplectic algorithms when numerically integrating Hamiltonian systems of differential equations. Runge--Kutta schemes are symplectic subject to a simple algebraic condition. Starting with Butcher's formalism it is shown that there exists a more natural basis for the set of necessary and sufficient order conditions for these methods, involving only s(s + 1)/2 free parameters for a symplectic s-stage scheme. A graph theoretical process for determining the new order conditions is outlined. Furthermore, it is shown that any rooted tree arising from the same free tree enforces the same algebraic constraint on the parametrized coefficients. When coupled with the standard simplifying assumptions for implicit schemes the number of order conditions may be further reduced. In the new framework a simple symmetry of the parameter matrix yields (not necessarily symplectic) self-adjoint methods. In this case the order conditions associated with even trees become redundant.

Order stars are a powerful modern tool for the development and analysis of numerical methods. They convey important information such as order and stability in a unified framework. A package for rendering order stars becomes part of the standard distribution in the next major release ofMathematica. An introduction to the theory is provided here, set in the context of numerical methods for Ordinary Differential Equations. The implementation is discussed and examples are given to illustrate why a computer algebra system is an ideal environment for the exploration of order stars.

Computer algebra provides a powerful tool for research in mathematics and the applied sciences. A symbolic package which may be used to explore and derive Runge-Kutta methods is presented. The package is the first of its kind to make use of the elegant graph-theoretical formalism attributed to Butcher.

It has been shown that numerical methods for Hamiltonian systems may be characterised in terms of so-called canonical elementary diierentials. Recent results by the authors demonstrate that the symplecticity condition for Runge-Kutta schemes may be eliminated by suitable parametrisation of the order conditions. These results are used to reexamine the relationship between order conditions for canonical Runge-Kutta methods and canonical elementary diierentials.