Variational system identification of the partial differential equations governing microstructure evolution in materials: Inference over sparse and spatially unrelated data

https://doi.org/10.1016/j.cma.2021.113706Get rights and content

Highlights

  • Our study focuses on inferring the equations that govern pattern formation.

  • We use a variational approach to system identification.

  • Our approach works with spatially unrelated, sparse and multi-source data.

  • We introduce a Confirmation Test of Consistency with Operator Suppression.

  • Synthetic data are used reflecting the characteristics of experimental images.

Abstract

Pattern formation is a widely observed phenomenon in diverse fields including materials physics, developmental biology and ecology, among many others. The physics underlying the patterns is specific to the mechanisms, and is encoded by partial differential equations (PDEs). With the aim of discovering hidden physics, we have previously presented a variational approach to identifying such systems of PDEs in the face of noisy data at varying fidelities (Computer Methods in Applied Mechanics and Engineering, 356:44–74, 2019). Here, we extend our variational system identification methods to address the challenges presented by image data on microstructures in materials physics. PDEs are formally posed as initial and boundary value problems over combinations of time intervals and spatial domains whose evolution is either fixed or can be tracked. However, the vast majority of microscopy techniques for evolving microstructure in a given material system deliver micrographs of pattern evolution over domains that bear no relation with each other at different time instants. The temporal resolution can rarely capture the fastest time scales that dominate the early dynamics, and noise abounds. Furthermore, data for evolution of the same phenomenon in a material system may well be obtained from different physical specimens. Against this backdrop of spatially unrelated, sparse and multi-source data, we exploit the variational framework to make judicious choices of weighting functions and identify PDE operators from the dynamics. A consistency condition arises for parsimonious inference of a minimal set of the spatial operators at steady state. It is complemented by a confirmation test that provides a sharp condition for acceptance of the inferred operators. The entire framework is demonstrated on synthetic data that reflect the characteristics of the experimental material microscopy images.

Introduction

Pattern-formation is ubiquitous in many branches of the physical sciences. It occurs prominently in material microstructures driven by diffusion, reaction and phase transformations, and is revealed by a range of microscopy techniques that delineate the components or constituent phases. In developmental biology, examples include the organization of cells in the early embryo, markings on animal coats, insect wings, plant petals and leaves, as well as the segregation of cell types during the establishment of tissues. In ecology, patterns are formed on larger scales as types of vegetation spreads across forests. For context, we briefly discuss the role of pattern forming systems of equations in these phenomena. Pattern formation during phase transformations in materials physics can happen as the result of instability-induced bifurcations from a uniform composition [1], [2], [3], which was the original setting of the Cahn–Hilliard treatment [4]. Following Alan Turing’s seminal work on reaction–diffusion systems [5], a robust literature has developed on the application of nonlinear versions of this class of partial differential equations (PDEs) to model pattern formation in developmental biology [6], [7], [8], [9], [10], [11], [12], [13], [14]. The Cahn–Hilliard phase field equation has also been applied to model other biological processes with evolving fronts, such as tumor growth and angiogenesis [15], [16], [17], [18], [19], [20], [21], [22]. Reaction–diffusion equations also appear in ecology, where they are more commonly referred to as activator–inhibitor systems, and are found to underlie large scale patterning [23], [24]. All of these pattern forming systems fall into the class of nonlinear, parabolic PDEs, and have spawned a vast literature in mathematical physics. They can be written as systems of first-order dynamics driven by a number of time-independent terms of algebraic and differential form. The spatio-temporal, differentio-algebraic operators act on either a composition (normalized concentration) or an order parameter. It also is common for the algebraic and differential terms to be coupled across multiple species.

Patterns in physical systems are of interest because, up to a point, human experts in each of the above fields (materials microscopists, developmental biologists and ecologists) are able to identify phenomena solely on the basis of patterns. This success of intuition fed by experience does, however, break down when, for instance, the materials scientist is confronted by the dynamic processes in an unstudied alloy, or the developmental biologist considers a previously neglected aspect of morphogenesis. In such settings the challenge is to discover the operative physics from among a range of mechanisms. As with all of quantitative physics, the only rigorous route to such discovery is the mathematical one. For systems that vary over space and time, this description is in the form of PDEs. Identification of participating PDE operators from spatio-temporal observations can thus uncover the underlying physical mechanisms, and lead to improved understanding, prediction and manipulation of these systems.

Distinct from classical adjoint-based approaches to inverse problems, the system identification problem on PDEs is to balance accuracy and complexity in finding the best model among many possible candidates that could explain the spatio-temporal data on the dynamical system. The proposed models can be parameterized by the coefficients of candidate operators, which serve as a basis. The task of solving this inverse problem can then be posed as finding the best coefficient values that allow a good agreement between model predictions and data. Sparsity of these coefficients is further motivated by the principle that parsimony of physical mechanisms is favored for applicability to wider regimes of initial and boundary value problems. The comparison procedure between models and data may be set up via the following two broad approaches. (a) When only sparse and noisy data are available, which might also be of indirect quantities that do not explicitly enter the PDEs, quantifying uncertainty in the identification becomes important and a Bayesian statistical approach is very useful. (b) When relatively larger volumes of observed data are accessible for quantities that directly participate in the PDEs, regression-based methods seeking to minimize an appropriate loss function can be highly efficient.

The first (Bayesian) approach has been successfully used for inferring a fixed number of coefficients in specified PDE models, leveraging efficient sampling algorithms such as Markov Chain Monte Carlo [25]. However, they generally require many forward solutions of the PDE models, which is computationally expensive and may quickly become impractical with growth in the number of PDE terms (i.e., dimension of the identification problem). As a result, the use of Bayesian approaches for identifying the best model from a large candidate set remains challenging. If made practical, they can be advantageous in providing uncertainty information, and offering significant flexibility in accommodating the data, which could be sparse and noisy, only available over some subset of the domain, collected at a few time steps, and composed of multiple statistical measures, functionals of the solution, or other indirect Quantities of Interests (QoIs). In the second (regression) approach, the PDEs themselves have to be represented by the data, which can be achieved by constructing the operators either in strong form such as finite difference, as in the Sparse Identification of Nonlinear Dynamics (SINDy) approach [26], or in weak form built upon basis functions, as in the Variational System Identification (VSI) approach [27]. Both representations pose stringent requirements on the data for accurately constructing the operators. The first of these is the need for time series data that can be related to chosen spatial points either directly or by interpolation. This is essential for consistency with a PDE that is written in terms of spatio-temporal operators at defined points in space and instants in time. The second requirement is for data with sufficient spatial resolution to construct spatial differential operators of possibly high order. However, these regression approaches enjoy the advantage over Bayesian methods that repeated forward PDE solutions are not necessary. PDE operators have thus been successfully identified from a comprehensive library of candidates [26], [27], [28], [29]. In a different approach to solving inverse problems [30], the strong form of a specified PDE is directly embedded in the loss function while training deep neural network representations of the solution variable. However, this approach depends on data at high spatial and temporal resolutions for successful training of the deep neural network representations of the solution variable. A perspective and comparison between these two approaches can be found in Ref. [31]. Lastly, there has been growing interest in combining graph theory with data-driven physics discovery. These approaches include genetic programming algorithm for discovering governing equations without pre-specification of the operators [32], and geometric learning for reconstructing the normal forms of equations using data collected from systems with unknown initial conditions and parameter values [33].

A significant discordance can exist in the form of material microstructure datasets when juxtaposed against the underlying premise of all the above approaches. These experimental datasets, while corresponding to different times are also commonly collected over different spatial subdomains of a physical specimen at each instant. This includes scenarios where the spatial subdomains have no overlap and are unrelated to each other. Furthermore, while subject to the same processing conditions, they may well come from different physical specimens. This is because the experimental techniques for extracting data at a given time after specimen preparation (which can include mechanical, chemical and thermal steps) involve destructive processes including cutting and grinding of the specimen. After this procedure, the entire specimen is removed from further experimentation, and a new specimen needs to be created if data is to be collected at a different time. Both of these conditions essentially negate the foundational notion of a PDE describing the temporal evolution of quantities at chosen spatial points for a single instantiation of that initial and boundary value problem (IBVP). The spatial subdomains on which microscopy is conducted may represent only small portions of entire specimens. Boundary data, if present, are also prey to the above loss of spatial localization over time. Finally, the effort of processing (specimen preparation, mechanical, chemical and heat treatment) to attain the desired kinetic rates and thermodynamic driving forces, and subsequently to obtain microscopy images, leads to sparsity of data in time; even tens of instants are uncommon.

In this communication, we extend our VSI methods [27] to sparse and spatially unrelated data, motivated by the above challenges of experimental microscopy in materials physics. The goal remains to discover the physical mechanisms underlying patterns of diffusion, reaction and phase transformation in materials physics. We present novel advances that exploit the variational framework and a notion of similarity between snapshots of data, which we make precise, to circumvent the loss of spatial relatedness so that temporally sparse data prove sufficient for inferring the kinetics. These methods are described in Section 2 and demonstrated on synthetic data in Section 2.6. Furthermore, imposition of consistency on the steady state versions of candidate PDEs opens the door to parsimonious choice of a minimal set of spatial operators. To complement this condition, we also present a confirmation test that is a sharp condition for acceptance of the inferred operators. The confirmation test is described in Section 3 and demonstrated on synthetic data in Section 3.2. Discussions and concluding remarks are presented in Section 4.

Section snippets

The Galerkin weak form for system identification

We first provide a brief discussion on the weak form of PDEs as used in this work. We start with the general strong form for first-order dynamics written as Ctχθ=0,where χ is a vector containing all possible linearly independent terms expressed as algebraic and differential operators on the scalar solution C: χ=[1,C,C2,,2C,],and θ is the vector of scalar coefficients of these operators. For example, using this nomenclature, the one-field diffusion reaction equation CtD2Cf=0with

Variational system identification with dynamic and steady state data

One challenge of the two-stage procedure introduced in Section 2 is data availability, where experimental measurements from the dynamic regime (non-steady state) of the system can be especially difficult and expensive to obtain. While the regression problems in Eqs. (40), (41) generally require over-determined systems to arrive at good solutions, this consequently demands a large dataset to delineate the relevant bases, since the number of rows in Eqs. (40), (41) is the number of time instants

Discussion and conclusions

The development of patterns in many physical phenomena is governed by a range of spatio-temporal PDEs. It is compelling to attempt to discover the analytic forms of these PDEs from data, because doing so immediately provides insight to the governing physics. System identification has been explored using the strong form [26], [28], [29], [43], [44] and the weak form [27], [34], [35], [36] of the PDEs as discussed in the Introduction. These techniques, however all need data that are spatially

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We acknowledge the support of Toyota Research Institute, United States of America, Award #849910, “Computational framework for data-driven, predictive, multi-scale and multi-physics modeling of battery materials” (ZW and KG). Additional support: This material is based upon work supported by the Defense Advanced Research Projects Agency (DARPA), United States of America under Agreement No. HR0011199002, “Artificial Intelligence guided multi-scale multi-physics framework for discovering complex

References (45)

  • DillonR. et al.

    Pattern formation in generalized turing systems i: Steady-state patterns in systems with mixed boundary conditions

    J. Math. Biol.

    (1994)
  • BarrioR.A. et al.

    A two-dimensional numerical study of spatial pattern formation in interacting turing systems

    Bull. Math. Biol.

    (1999)
  • BarrioR.A. et al.

    Modeling the skin pattern of fishes

    Phys. Rev. E

    (2009)
  • MainiPhilip K. et al.

    Turing’s model for biological pattern formation and the robustness problem

    Interface Focus

    (2012)
  • SpillF. et al.

    Hybrid approaches for multiple-species stochastic reaction–diffusion models

    J. Comput. Phys.

    (2015)
  • KorvasováK. et al.

    Investigating the turing conditions for diffusion-driven instability in the presence of a binding immobile substrate

    J. Theoret. Biol.

    (2015)
  • GarikipatiK.

    Perspectives on the mathematics of biological patterning and morphogenesis

    J. Mech. Phys. Solids

    (2017)
  • WiseS.M. et al.

    Three-dimensional multispecies nonlinear tumor growth–model and numerical method

    J. Theoret. Biol.

    (2008)
  • CristiniV. et al.

    Nonlinear simulations of solid tumor growth using a mixture model: invasion and branching

    J. Math. Biol.

    (2009)
  • LowengrubJ.S. et al.

    Nonlinear modelling of cancer: bridging the gap between cells and tumours

    Nonlinearity

    (2010)
  • LowengrubJ.S. et al.

    Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission

    Phys. Rev. E

    (2009)
  • VilanovaG. et al.

    Capillary networks in tumor angiogenesis: From discrete endothelial cells to phase-field averaged descriptions via isogeometric analysis

    Numer. Methods Biomed. Eng.

    (2013)
  • Cited by (0)

    View full text