Paper deep dive
Inelastic Constitutive Kolmogorov-Arnold Networks: A generalized framework for automated discovery of interpretable inelastic material models
Chenyi Ji, Kian P. Abdolazizi, Hagen Holthusen, Christian J. Cyron, Kevin Linka
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 7/20/2026, 11:49:52 PM
Summary
The paper introduces inelastic Constitutive Kolmogorov-Arnold Networks (iCKANs), a novel neural network architecture that automates the discovery of interpretable, closed-form symbolic constitutive laws for inelastic materials. By integrating Kolmogorov-Arnold Networks with generalized inelastic constitutive frameworks, iCKANs can extract elastic and inelastic potential functions from experimental data (e.g., stress-strain curves) while preserving physical interpretability and handling additional features like temperature.
Entities (10)
Relation Signals (9)
iCKANs → discovers → symbolic constitutive laws
confidence 95% · This novel artificial neural network architecture can discover in an automated manner symbolic constitutive laws describing both the elastic and inelastic behavior of materials.
Chenyi Jia → affiliatedwith → RWTH Aachen University
confidence 90% · Chenyi Jia... Computational Mechanics in Medicine, Applied Medical Engineering, RWTH Aachen University
Kian P. Abdolazizi → affiliatedwith → Hamburg University of Technology
confidence 90% · Kian P. Abdolazizi... Institute for Continuum and Material Mechanics, Hamburg University of Technology
iCKANs → extends → CKANs
confidence 90% · In this work, we advance KAN-based constitutive modeling by extending it to inelastic materials and introducing inelastic Constitutive Kolmogorov–Arnold Networks (iCKANs).
iCKANs → outputs → Elastic Potential
confidence 90% · That is, it can translate data from material testing into corresponding elastic and inelastic potential functions in closed mathematical form.
iCKANs → outputs → Inelastic Potential
confidence 90% · That is, it can translate data from material testing into corresponding elastic and inelastic potential functions in closed mathematical form.
iCKANs → processes → VHB 4905
confidence 90% · We demonstrate the advantages of iCKANs using both synthetic data and experimental data of the viscoelastic polymer materials VHB 4910 and VHB 4905.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:A key problem of solid mechanics is the identification of the constitutive law of a material, that is, the relation between strain history and stress. Machine learning has lead to considerable advances in this field lately. Here we introduce inelastic Constitutive Kolmogorov-Arnold Networks (iCKANs). This novel artificial neural network architecture can discover in an automated manner symbolic constitutive laws describing both the elastic and inelastic behavior of materials. That is, it can translate data from material testing into corresponding elastic and inelastic potential functions in closed mathematical form. We demonstrate the advantages of iCKANs using both synthetic data and experimental data of the viscoelastic polymer materials VHB 4910 and VHB 4905. The results demonstrate that iCKANs accurately capture complex viscoelastic behavior while preserving physical interpretability. It is a particular strength of iCKANs that they can process not only mechanical data but also arbitrary additional information available about a material (e.g., about temperature-dependent behavior). This makes iCKANs a powerful tool to discover in the future also how specific processing or service conditions affect the properties of materials.
Tags
Links
- Source: https://arxiv.org/abs/2602.17750v3
- Canonical: https://arxiv.org/abs/2602.17750v3
Trouble viewing inline? Open PDF directly →
Full Text
131,333 characters extracted from source content.
Expand or collapse full text
Inelastic Constitutive Kolmogorov-Arnold Networks: A generalized framework for automated discovery of interpretable inelastic material models Chenyi JiaChenyi Ji^a, Kian P. AbdolazizibKian P. Abdolazizi^b, Hagen HolthusencHagen Holthusen^c, Christian J. Cyronb,dChristian J. Cyron^b,d, Kevin Linkaa,∗ a Computational Mechanics in Medicine, Applied Medical Engineering, RWTH Aachen University Pauwelsstraße 20, 52074 Aachen, Germany b Institute for Continuum and Material Mechanics, Hamburg University of Technology Eißendorfer Straße 42, 21073 Hamburg, Germany c Institute of Applied Mechanics, University of Erlangen-Nuremberg Egerlandstrasse 5, 91058 Erlangen, Germany d Institute of Material Systems Modeling, Helmholtz-Zentrum Hereon Max-Planck-Straße 1, 21502 Geesthacht, Germany Abstract. A key problem of solid mechanics is the identification of the constitutive law of a material, that is, the relation between strain and stress. Machine learning has lead to considerable advances in this field lately. Here we introduce inelastic Constitutive Kolmogorov–Arnold Networks (iCKANs). This novel artificial neural network architecture can discover in an automated manner symbolic constitutive laws describing both the elastic and inelastic behavior of materials. That is, it can translate data from material testing into corresponding elastic and inelastic potential functions in closed mathematical form. We demonstrate the advantages of iCKANs using both synthetic data and experimental data of the viscoelastic polymer materials VHB 4910 and VHB 4905. The results demonstrate that iCKANs accurately capture complex viscoelastic behavior while preserving physical interpretability. It is a particular strength of iCKANs that they can process not only mechanical data but also arbitrary additional information available about a material (e.g., about temperature-dependent behavior). This makes iCKANs a powerful tool to discover in the future also how specific processing or service conditions affect the properties of materials. 11footnotetext: Corresponding author linka@ame.rwth-aachen.de Keywords: Kolmogorov-Arnold Networks, constitutive modeling, inelasticity, finite strains, model discovery, symbolic regression 1 Introduction Accurately predicting the behavior of complex material systems is a core challenge in engineering and the physical sciences. Conventional computational mechanics relies on constitutive models to describe how materials respond to external stimuli. However, traditional constitutive models, grounded in continuum mechanics, depend on simplifying assumptions that often fail to capture the full complexity of real-world material behavior—particularly for inelastic materials. Moreover, developing these models is typically incremental, labor-intensive, and demands substantial domain expertise, with the process repeated for each novel material system. To overcome these limitations, data-driven material modeling has emerged as a promising alternative, bypassing explicit constitutive model formulation and leveraging experimental data directly [Fuhg et al., 2024]. The model-free or direct data-driven paradigm, pioneered by Kirchdoerfer and Ortiz [2016], eliminates predefined constitutive equations in favor of using raw material measurements. Another avenue is symbolic regression, which constructs closed-form mathematical models from data by systematically combining analytical expressions to achieve a balance between accuracy and interpretability [Versino et al., 2017; Bomarito et al., 2021; Kabliman et al., 2021; Abdusalamov et al., 2023]. Expanding on this, EUCLID (Efficient Unsupervised Constitutive Law Identification and Discovery) applies sparse regression to a broad model library, enabling learning from full-field data [Flaschel et al., 2021; Flaschel, 2023; Flaschel et al., 2023]. Recent developments have further enhanced these approaches to better handle noisy experimental data [Narouie et al., 2026]. A distinct class of data-driven approaches involves constitutive neural networks, which have gained prominence due to their flexibility and capacity for universal function approximation [Hornik et al., 1989]. Unlike methods that impose physical constraints via the loss function [Raissi et al., 2019; Masi et al., 2021], these networks encode physical consistency directly through architectural design. Notable variants include Constitutive Artificial Neural Networks (CANNs) [Linka et al., 2021; Linka and Kuhl, 2023], which favor interpretable, sparse architectures, and Physics-Augmented Neural Networks (PANNs) [Rosenkranz et al., 2024; Linden et al., 2023], which leverage denser, more expressive structures. Additional neural network classes have further broadened the modeling landscape: neural ordinary differential equations (ODEs) have been used to identify polyconvex strain energy functions [Tac et al., 2022]; graph neural networks [Maurizi et al., 2022] and neural operators [You et al., 2022] have been explored for learning complex constitutive and surrogate models. These methods have also found applications in the design and optimization of metamaterials [Fernández et al., 2021, 2022]. Current research pushes the boundaries by targeting increasingly complex phenomena, such as anisotropy and inelasticity [Holthusen and Kuhl, 2026]. To capture inelastic effects, the Generalized Standard Materials framework offers a broad and physically grounded foundation [Halphen and Nguyen, 1975]. This approach introduces an inelastic dissipation potential alongside the elastic free energy, governing the evolution of internal variables and inelastic deformation in a thermodynamically consistent manner. Integrating this framework with constitutive neural networks has enabled the discovery of material laws for finite-strain viscoelasticity [Rosenkranz et al., 2024; As’ ad and Farhat, 2023; Holthusen et al., 2026; Taç et al., 2023], plasticity [Flaschel et al., 2025; Boes et al., 2026; Jadoon et al., 2025], fracture [Dammaß et al., 2025], and biological processes like growth and remodeling [Holthusen et al., 2025]. Constitutive models based on internal variables have also been successfully employed in neural ODEs and neural operators to model history- and path-dependent behaviors [Jones et al., 2022; Guo et al., 2025]. Furthermore, approaches such as long short-term memory networks have been introduced to enhance computational efficiency in plasticity modeling [Li et al., 2025]. Constitutive Kolmogorov–Arnold Networks (CKANs) [Abdolazizi et al., 2025; Thakolkaran et al., 2025] offer a distinctive alternative by employing B-Splines as nonlinear, trainable activation functions whose mathematical forms can be extracted for interpretation [Liu et al., 2024]. The introduction of Kolmogorov–Arnold Networks (KANs) has prompted comparative studies with traditional neural networks across a range of applications. KANs have shown improved accuracy and convergence in solving partial differential equations relative to multilayer perceptrons [Wang et al., 2025; Abueidda et al., 2025; Kiyani et al., 2025], with their greatest strengths evident in symbolic regression tasks [Ji et al., 2024]. In the context of material modeling, KANs have been incorporated into CKANs and Input-Convex KANs to enhance interpretability and facilitate the discovery of convex strain energy functions, all while retaining the adaptability of neural network-based frameworks [Abdolazizi et al., 2025; Thakolkaran et al., 2025]. These developments underscore the potential of KAN-based architectures to increase the transparency of data-driven material models without sacrificing predictive accuracy, though challenges in modeling inelastic behavior persist. In this work, we advance KAN-based constitutive modeling by extending it to inelastic materials and introducing inelastic Constitutive Kolmogorov–Arnold Networks (iCKANs). This novel approach integrates the generalized inelastic constitutive framework with the inherent interpretability of KANs, aiming to deliver a robust, data-driven, and transparent model for complex inelastic behavior. iCKANs leverage experimental stress–strain data, optionally enriched with non-mechanical features, to uncover the underlying elastic and inelastic potentials of materials. Crucially, the trainable activations of the KANs are subsequently symbolized, yielding interpretable closed-form expressions for these potentials. By enabling the automated derivation of symbolic inelastic potentials, previously accessible primarily to a limited group of experts, iCKANs address a longstanding challenge in finite strain inelasticity. The overall workflow of the iCKAN framework is illustrated in Figure 1. Outline. Section 2 provides a concise review of the generalized constitutive framework for inelastic materials at finite strains. Section 3 details the proposed iCKANs’ methodology, beginning with the foundations of Kolmogorov–Arnold Networks and their partially input-convex variants, followed by their application to constitutive modeling. Section 4 addresses the symbolic extraction process to obtain interpretable, closed-form representations of the discovered models. Section 5 evaluates the performance of iCKANs on both synthetic and experimental viscoelastic datasets. The paper concludes in Section 6 with a discussion of key findings and future research directions. Experimental dataKolmogorov-Arnold NetworksSymbolic function ∙ Time-dependent stress- deformation data ∙ Non-mechanical data f =1f= [rgb]0.96484375,0.66015625,0 [named]pgfstrokecolorrgb0.96484375,0.66015625,0f_1=2f= [rgb]0,0.328125,0.625 [named]pgfstrokecolorrgb0,0.328125,0.625f_2=3f= [rgb]0.33984375,0.671875,0.15234375 [named]pgfstrokecolorrgb0.33984375,0.671875,0.15234375f_3Stretch FiF_i Stress PiP_i Training with sparsification Stress PDeformation gradient FElastic pot.Inelastic pot. ∙ Recurrent network architecture ∙ Mathematical/physical constraints ∙ (Partially) convex splines Symbolification Stress predictionsP=a(F+b)2P=a(F+b)^2=1f= [rgb]0.96484375,0.66015625,0 [named]pgfstrokecolorrgb0.96484375,0.66015625,0f_1=2f= [rgb]0,0.328125,0.625 [named]pgfstrokecolorrgb0,0.328125,0.625f_2=3f= [rgb]0.33984375,0.671875,0.15234375 [named]pgfstrokecolorrgb0.33984375,0.671875,0.15234375f_3Stretch FiF_i Stress PiP_i ∙ Interpretability of elastic and inelastic potential ∙ Efficient numerical simulations Figure 1: Inelastic Constitutive Kolmogorov-Arnold Network (iCKAN) pipeline for automated interpretable model discovery of inelastic materials. Three-dimensional stress-strain data with optional additional features (e.g., temperature) collected in a feature vector f are used to train the iCKAN model, which consists of two KANs representing the elastic potential ψ and the inelastic potential ω. The trained model can then be analyzed using symbolic regression to extract interpretable mathematical expressions. 2 Constitutive modeling of inelastic materials In this section, we briefly review the generalized constitutive framework for inelastic materials at finite strains by Holthusen et al. [2023, 2026], which can be seen as an equivalence to the framework of Generalized Standard Materials [Halphen and Nguyen, 1975]. The framework employs the multiplicative decomposition of the deformation gradient and is based on two fundamental scalar-valued quantities, the elastic potential ψ and an inelastic potential ω [Holthusen et al., 2024b]. Multiplicative decomposition. The general assumption of multiplicative decomposition is applied on the deformation gradient F, resulting into an elastic part eF_e and inelastic part iF_i, i.e., =eiF=F_eF_i. Conceptually, an intermediate configuration is introduced, which is reached by the inelastic deformation iF_i from the reference configuration. The elastic deformation eF_e then maps the intermediate configuration to the current configuration. However, the intermediate configuration is fictitious and not unique, since any rotation †∈SO(3)Q ∈ SO(3) applied to the intermediate configuration does not alter the physics of the resulting deformation gradient, i.e., =e†T†i=e†i†F=F_eQ ^TQ F_i=F_e F_i . Following Holthusen et al. [2023, 2024b], a unique co-rotated intermediate configuration is obtained by applying the rotation of the polar decomposition of i=iiF_i=R_iU_i to the intermediate configuration, i.e., ¯e=ei=i−1 F_e=F_eR_i=FU_i^-1, with iU_i and iR_i being the stretch part and rotation part of the inelastic deformation, respectively. Fi†F_i e†F_e iF_ieF_e†Q iU_iiR_ireferenceconfigurationintermediate configurationcurrent configurationarbitrarily rotatedintermediate configurationco-rotatedintermediate configurationℬ0B_0ℬtB_t Figure 2: Multiplicative split of the deformation gradient into elastic and inelastic part, including the non-uniqueness of the intermediate configuration and the co-rotated intermediate configuration [Holthusen et al., 2023]. The right Cauchy-Green deformation tensors are defined as =TC=F^TF and ¯e=i−1i−1 C_e=U_i^-1CU_i^-1 in the current and co-rotated intermediate configuration, respectively. Due to the fact that ¯e=i−1i−1 C_e=U_i^-1CU_i^-1 serves as an unique measurement for elastic stretches, the elastic potential can be defined depending on it, i.e., ψ=ψ(¯e)ψ=ψ( C_e), while preserving the principles of objectivity and frame indifference. To ensure that the elastic potential is a scalar-valued isotropic function, we define it based on the principal invariants of ¯e C_e, i.e., ψ=ψ(I1¯e,I2¯e,I3¯e)ψ=ψ(I_1 C_e,I_2 C_e,I_3 C_e), with I1¯e=tr(¯e),I2¯e=12[(tr(¯e))2−tr(¯e2)],I3¯e=det(¯e).I_1 C_e=tr( C_e), I_2 C_e= 12 [(tr( C_e))^2-tr( C_e^2) ], I_3 C_e= ( C_e)\,. (1) Dissipation inequality. Following the second law of thermodynamics, the dissipation D has to be non-negative so that thermodynamic consistency is preserved. This can be expressed using Clausius-Plank inequality =:12˙−ψ˙≥0D=S: 12 C- ψ≥ 0. By inserting the fact that the elastic potential depends on ¯e C_e and applying the chain rule, we obtain =(−2i−1∂ψ∂¯ei−1):12˙+2¯e∂ψ∂¯e:˙ii−1≥0,D= (S-2\,U_i^-1 ∂ψ∂ C_eU_i^-1 ): 12 C+2\, C_e ∂ψ∂ C_e: U_iU_i^-1≥ 0\,, (2) which has to be fulfilled for any arbitrary value of ˙ C. Hence, the second Piola-Kirchhoff stress tensor and the symmetric elastic Mandel-like stress tensor are defined as =2i−1∂ψ∂¯ei−1and¯=2¯e∂ψ∂¯e,S=2\,U_i^-1 ∂ψ∂ C_eU_i^-1 =2\, C_e ∂ψ∂ C_e\,, (3) respectively. The first Piola-Kirchhoff stress can be derived by pulling back the second Piola-Kirchhoff stress, i.e., =P=FS. Moreover, by noting that ¯ is symmetric, the dissipation inequality can be reduced to =¯:˙ii−1=¯:¯i≥0with¯i=sym(˙ii−1).D= : U_iU_i^-1= : D_i≥ 0 D_i=sym( U_iU_i^-1)\,. (4) Evolution equation. At this stage, a suitable evolution equation for the inelastic rate tensor ¯i D_i remains to be specified. To this end, a dual inelastic potential ω is introduced, which is assumed to be a scalar-valued, isotropic function of the thermodynamically consistent driving force ¯ , i.e., ω=ω(¯)ω=ω( ). It has been shown that if ω is constructed to be convex, zero-valued at ¯= =0, and non-negative, the reduced dissipation inequality is satisfied automatically [Germain et al., 1983]. While these conditions are sufficient, they are not necessary to guarantee non-negative dissipation. Following [Holthusen et al., 2026], a more general formulation is obtained by introducing an invariant-based inelastic potential of the form ω(¯)=ω(I1¯,J2¯,J3¯3),ω( )=ω\! (I_1 , J_2 , [3]J_3 ), (5) which is required to be convex, zero-valued at its origin, and non-negative with respect to the three invariants***It is worth noting that the third invariant is generally not convex with respect to ¯ . I1¯ I_1 =tr(¯), =tr( ), J2¯ J_2 =12tr(dev(¯)2), = 12tr\! (dev( )^2 ), J3¯ J_3 =13tr(dev(¯)3), = 13tr\! (dev( )^3 ), (6) ∂I1¯∂¯ ∂ I_1 ∂ =, =I, ∂J2¯∂¯ ∂ J_2 ∂ =dev(¯), =dev( ), ∂J3¯∂¯ ∂ J_3 ∂ =dev(dev(¯)2). =dev\! (dev( )^2 ). Here, dev(¯)=¯−13I1¯dev( )= - 13I_1 I denotes the deviatoric part of a second-order tensor. It is worth emphasizing that the square and cubic roots of the second and third invariants, respectively, are essential to ensure non-negative dissipation within this framework [Holthusen et al., 2026]. The evolution equation for the inelastic rate tensor then follows directly from the inelastic potential as ¯i=∂ω(¯)∂¯=∂ω∂I1¯+∂ω∂J2¯dev(¯)+∂ω∂J3¯dev(dev(¯)2). D_i= ∂ω( )∂ = ∂ω∂ I_1 \,I+ ∂ω∂ J_2 \,dev( )+ ∂ω∂ J_3 \,dev\! (dev( )^2 ). (7) Incompressible materials. For the special case of incompressible materials, i.e., det()=1 (F)=1, a Lagrange multiplier term is added to the elastic potential ψ to obtain the augmented elastic potential ψ ψ, ψ^=ψ(¯e)−p(I3). ψ=ψ( C_e)-p(I_3^C)\,. (8) The term p can be considered as the hydrostatic pressure and is determined using the boundary conditions. It is important to note that the elastic part of the deformation is, in general, not incompressible, meaning that I3¯eI_3 C_e not necessarily remain equal to one [Holthusen et al., 2024b]. 3 Inelastic Constitutive Kolmogorov–Arnold Networks (iCKANs) Building on the constitutive framework presented in the previous section, this section extends the idea of CKANs [Abdolazizi et al., 2025] by introducing inelastic Constitutive Kolmogorov–Arnold Networks (iCKANs) for automatic model discovery of inelastic materials. We start with a brief review of Kolmogorov-Arnold Networks (KANs) [Liu et al., 2024] and their modified variant of monotonic input-convex KANs [Thakolkaran et al., 2025], and then we continue with presenting partially input-convex KANs, in which the output is generally convex with respect to specified inputs. Afterwards, we apply these partially input-convex KANs to express the elastic and inelastic potential functions of the constitutive material model introduced in the last section, leading to the formulation of iCKANs. Thus, iCKANs naturally preserve mandatory constraints and flexibility due to the expression of activation functions via B-Splines. 3.1 Input-convex Kolmogorov-Arnold Networks Serving as a promising alternative to mutli-layer perceptrons, KANs replace the combination of a linear layer following with a nonlinear activation function with a naturally nonlinear, trainable activation. KANs are inspired by the Kolmogorov-Arnold representation theorem, which states that any continuous multivariate function on a bounded domain f:[0,1]n→ℝf:[0,1]^n can be represented as a finite composition of univariate functions f(0)=f(x0,1,x0,2,…,x0,n)=∑j=12n+1ϕ1,1,j(∑i=1nϕ0,j,i(x0,i)),f(x_0)=f(x_0,1,x_0,2,…,x_0,n)= _j=1^2n+1 _1,1,j ( _i=1^n _0,j,i(x_0,i) )\,, (9) where x0,ix_0,i are the elements of the input vector 0x_0 and ϕ0,j,i:[0,1]→ℝ _0,j,i:[0,1] and ϕ1,1,j:ℝ→ℝ _1,1,j:R are univariate continuous functions [Kolmogorov, 1961; Abdolazizi et al., 2025]. This equation can be seen as a two-layer network with a topology of [n,2n+1,1][n,2n+1,1]. Extending this idea, a deeper architecture with L layers can be constructed, resulting in a KAN with the topolgy of [n0,n1,…,nL][n_0,n_1,…,n_L], where nin_i is the number of neurons in the i-th layer. Thus, KANs can be expressed as the following nested summation f()=∑iL−1=1nL−1ϕL−1,iL,iL−1(∑iL−2=1nL−2…(∑i2=1n2ϕ2,i3,i2(∑i1=1n1ϕ1,i2,i1(∑i0=1n0ϕ0,i1,i0(xi0))))…).f(x)= _i_L-1=1^n_L-1 _L-1,i_L,i_L-1 ( _i_L-2=1^n_L-2… ( _i_2=1^n_2 _2,i_3,i_2 ( _i_1=1^n_1 _1,i_2,i_1 ( _i_0=1^n_0 _0,i_1,i_0(x_i_0) ) ) )… )\,. (10) The activation function connecting the i-th neuron in the l-th layer (l,i)(l,i) and the i-th neuron in the (l+1)(l+1)-th layer (l+1,j)(l+1,j) is denoted as ϕl,j,i _l,j,i for l=0,…,L−1l=0,…,L-1, i=1,…,nli=1,…,n_l and j=1,…,nl+1j=1,…,n_l+1. An activation function is defined as follows ϕl,j,i(x)=wb⋅b(x)+ws⋅s(x) _l,j,i(x)=w_b· b(x)+w_s· s(x) (11) with a basis function b(x)b(x), usually b(x)=x/(1+e−x)b(x)=x/(1+e^-x), and a spline function s(x)=∑iciBi(x)s(x)= _ic_iB_i(x), where cic_i are trainable control points and BiB_i are B-spline basis functions. The factors wbw_b and wsw_s are trainable coefficients. The output of neuron (l+1,j)(l+1,j), xl+1,jx_l+1,j, is defined as the sum of the inputs xl,ix_l,i after transformation, i.e., xl+1,j=∑i=1nlϕl,j,i(xl,i)x_l+1,j= _i=1^n_l _l,j,i(x_l,i), where each input xl,ix_l,i is processed by its associated univariate activation function ϕl,j,i _l,j,i. For the simplicity of notation, the layer-wise mapping is defined using the nonlinear function matrix l _l as follows l+1=l(l)=[ϕl,1,1(∙)ϕl,1,2(∙)…ϕl,1,nl(∙)ϕl,2,1(∙)ϕl,2,2(∙)…ϕl,2,nl(∙)⋮⋱⋮ϕl,nl+1,1(∙)ϕl,nl+1,2(∙)…ϕl,nl+1,nl(∙)]l.x_l+1= _l(x_l)= bmatrix _l,1,1( )& _l,1,2( )&…& _l,1,n_l( )\\ _l,2,1( )& _l,2,2( )&…& _l,2,n_l( )\\ & & & \\ _l,n_l+1,1( )& _l,n_l+1,2( )&…& _l,n_l+1,n_l( ) bmatrixx_l\,. (12) Thus, Equation 10 can be written compactly as KAN()=(ΦL−1∘ΦL−2∘⋯∘Φ1∘Φ0)().KAN(x)=( _L-1 _L-2 … _1 _0)(x)\,. (13) 3.1.1 Fully input-convex architecture A function that is twice differentiable is convex if and only if its second derivative is nonnegative over its entire domain. For instance, consider a two-layer KAN with the input 0=x0,1,…,x0,n0x_0=\x_0,1,…,x_0,n_0\, we can express it as the following function f(0)=∑j=1niϕ1,1,j(∑i=1n0ϕ0,j,i(x0,i)).f(x_0)= _j=1^n_i _1,1,j ( _i=1^n_0 _0,j,i(x_0,i) )\,. (14) The first and second derivatives of f with respect to x0,kx_0,k are given by ∂f∂x0,k=∑j=1n1ϕ1,1,j′(∑i=1n0ϕ0,j,i(x0,k))⋅ϕ0,j,k′(x0,k), ∂ f∂ x_0,k= _j=1^n_1 _1,1,j ( _i=1^n_0 _0,j,i(x_0,k) )· _0,j,k (x_0,k)\,, (15) ∂2f∂x0,k2=∑j=1n1ϕ1,1,j′(∑i=1n0ϕ0,j,i(x0,i))⋅ϕ0,j,k′(x0,k)2+∑j=1n1ϕ1,1,j′(∑i=1n0ϕ0,j,i(x0,i))⋅ϕ0,j,k′(x0,k), ∂^2f∂ x_0,k^2= _j=1^n_1 _1,1,j ( _i=1^n_0 _0,j,i(x_0,i) )· _0,j,k (x_0,k)^2+ _j=1^n_1 _1,1,j ( _i=1^n_0 _0,j,i(x_0,i) )· _0,j,k (x_0,k)\,, (16) respectively. Consequently, a sufficient condition for f to be convex with respect to the inputs 0x_0 is that the following criteria hold ϕ0,j,i′≥0∀i=1,…,n0;∀j=1,…,n1;ϕ1,1,j′≥0∀j=1,…,n1;ϕ1,1,j′≥0∀j=1,…,n1. split _0,j,i &≥ 0 ∀ i=1,…,n_0\,;\,∀ j=1,…,n_1\,;\\ _1,1,j &≥ 0 ∀ j=1,…,n_1\,;\\ _1,1,j &≥ 0 ∀ j=1,…,n_1\,. split (17) That is, all first-layer activation functions ϕ0,j,i _0,j,i must be convex, while all second-layer activation functions ϕ1,1,j _1,1,j must be both convex and non-decreasing. Convex and non-decreasing activations. Constraining all activations to be convex and monotonic non-decreasing results into a monotonic input-convex KAN. To satisfy the aforementioned criteria and achieve monotonically increasing and convex activations, Thakolkaran et al. [2025] proposes to constrain the control points of the B-splines activation functions. First, zero base functions are chosen, i.e. b(x)=0b(x)=0. Then, wsw_s is constrained to be positive, monotonically increasing, and convex by ws∗=softplus(ws)w_s^*=softplus(w_s). Finally, the control points cic_i are modified to become convex coefficients ci∗c_i^*. Thus, s∗(x)=∑ici∗Bi(x)s^*(x)= _ic_i^*B_i(x) is monotonically increasing and convex, i.e., ci+2−ci+1≥ci+1−ci≥0,∀i.c_i+2-c_i+1≥ c_i+1-c_i≥ 0,\,∀ i\,. (18) Finally, the original activation in Equation 11 is modified to ϕ(x)=ws∗⋅∑ici∗Bi(x).φ(x)=w_s^*· _ic_i^*B_i(x)\,. (19) To ensure that the activations outside the initial grid range also remain convex and monotonically increasing, linear extrapolation is applied at both ends of the grid range [Polo-Molina et al., 2024; Thakolkaran et al., 2025]. Generally convex activations. The previously introduced approach provides a convex and non-decreasing KAN-output with respect to its inputs. However, in some cases, specifically for the inelastic potential, the output is required to be convex with a desired stationary point. Therefore, we apply an additional activation function to the KAN-output to ensure the convexity while allowing for non-monotonic behavior. It is important to note that the network architecture remains unchanged from the monotonic input-convex KANs and only the output of the KAN is passed through this additional activation function. Given a convex and non-decreasing multivariate function f()f(x), one can construct a new function f~() f(x) that is convex, zero-valued at its origin and has a zero derivative at =αx=α by passing it through the following defined operation ℋ(∙)H( ), f~()=ℋ(f())=f()−f(α)−∇f()|=α⋅(−α), f(x)=H(f(x))=f(x)-f(α)- .\ (x) |_x=α·(x-α)\,, (20) For a more general case, we can apply a further manipulation step, f^()=max(f~()−c,0)with c>0. f(x)=max( f(x)-c,0) c>0\,. (21) Applying the resulting functions f~() f(x) and f^() f(x) to the KAN-output preserve convexity. A demonstration of this for an univariate simplification is shown in 3(a). The proof of this theorem can be found in Appendix A.1. 3.1.2 Partially input-convex architecture In general, the outputs are required to be convex with respect to selected input arguments, but not necessarily with respect to all inputs. In particular, the elastic potential and the inelastic potential must be convex in the invariant-based arguments, while convexity with respect to additional non-mechanical features is not required. For this reason, we developed a partially input-convex KAN, as shown in 3(b) [Amos et al., 2017; Deschatre and Warin, 2025]. The output z=f(,)z=f(x,y) is convex with respect to the input variables =[x1,x2,x3]x=[x_1,x_2,x_3], but not necessarily convex with respect to the input variable =[y1]y=[y_1]. This is achieved by modifying the first layer of a generally input-convex KAN. While the activations corresponding to the convex input variables x remain convex and non-decreasing, the activations with respect to the input variable y are unconstrained. The outputs of the first layer are then additively combined and passed through the second layer. It is to note, that this implementation is a simplified variant of partially input-convex networks [Amos et al., 2017], as a more generalized implementation is beyond the scope of this work. (a) Postprocessing of a convex, non-decreasing function x1x_1x2x_2x3x_3y1y_1z (b) A partially input-convex Kolmogorov-Arnold Network Figure 3: (a) Demonstration of the postprocessing of a convex and non-decreasing function f(x)f(x) to achieve a general convex function with the designed stationary point at x=0x=0. f~(x) f(x) is achieved by ℋH-operation in Equation 20 and f^(x) f(x) is achieved by Equation 21. (b) A partially input-convex Kolmogorov-Arnold Network, where the output z is convex with respect to the yellow-marked input variables (x1,x2,x3)(x_1,x_2,x_3) and unconstrained to the blue-marked input variables y1y_1. 3.2 Constitutive modeling with KANs Recalling the constitutive framework in Section 2, the elastic potential ψ and inelastic potential ω must be defined to formulate an inelastic material. In this work, we propose using partially input-convex KANs to model these two scalar-valued functions. Thus, we present the inelastic Constitutive Kolmogorov-Arnold Networks (iCKANs) framework, which can be considered an extension to the Constitutive Kolmogorov-Arnold Networks (CKANs) for hyperelastic materials introduced in Abdolazizi et al. [2025]. 3.2.1 Elastic potential CKANs have been formulated using different choices of functional bases for the elastic potential function, e.g., pricipal invariants and principal stretches [Abdolazizi et al., 2025]. In this work, we adopt the variant, in which the principal invariants are selected as the functional basis for constructing the KAN representing the elastic potential. In the original CKAN implementation, an input-monotonic KAN is employed [Polo-Molina et al., 2024]. Here, we use the previously introduced input-convex KAN to enforce convexity of the learned elastic potential representation. The elastic potential function ψ should be constructed polyconvex with respect to the deformation gradient F, to ensure the existence of minimizers of the elastic potential function. To achieve a polyconvex elastic potential function, the input arguments of the KAN must be polyconvex in F [Hartmann and Neff, 2003]. It is evident that the elastic potential should be dependent on the principal invariants of the elastic right Cauchy-Green deformation tensor in the co-rotated intermediate configuration ¯e=i−1i−1 C_e=U_i^-1CU_i^-1. For nearly incompressible materials, it is common practice to decouple the volumetric and isochoric parts of the elastic deformation [Holthusen et al., 2024b; Flory, 1961]. Therefore, the polyconvex modified elastic invariants (I^1¯e,I^2¯e,I^3¯e)( I_1 C_e, I_2 C_e, I_3 C_e) are then defined using the principal invariants (Equation 1) as [Hartmann and Neff, 2003; Holthusen et al., 2024a] I^1¯e=I~1¯e−3,I^2¯e=(I~2¯e)3/2−33/2,I^3¯e=(J¯e−1)2, I_1 C_e= I_1 C_e-3, I_2 C_e=( I_2 C_e)^3/2-3^3/2, I_3 C_e=(J C_e-1)^2\,, (22) withJ¯e=det(¯e),I~1¯e=(J¯e)−2/3I1¯e,I~2¯e=(J¯e)−4/3I2¯e. J C_e= ( C_e)\,, I_1 C_e= (J C_e )^-2/3I_1 C_e, I_2 C_e= (J C_e )^-4/3I_2 C_e\,. (23) Thus, a polyconvex elastic potential is constructed using the modified invariants of ¯e C_e as input argument for the input-convex KAN ψ=ψ(¯e)=ψ(I^1¯e,I^2¯e,I^3¯e).ψ=ψ( C_e)=ψ( I_1 C_e, I_2 C_e, I_3 C_e)\,. (24) It is to note that the polyconvexity with respect to the arguments is not a necessary condition for the elastic potential to fulfill the physics, but rather a convenient assumption for the potential minimization. Furthermore, a proof that the output elastic potential ψ and its derivative are zero at undeformed state can be found in Appendix A.2 3.2.2 Inelastic potential What remains to be defined is the inelastic potential ω, which governs the evolution equation of the inelastic strain. To satisfy thermodynamic consistency, the inelastic potential ω must be zero-valued at its origin, nonnegative and convex with respect to its arguments, which are the modified stress invariants of the elastic Mandel stress ¯ (Equation 6), I^1¯=I1¯,J^2¯=J2¯,J^3¯=J3¯3. I_1 =I_1 , J_2 = J_2 , J_3 = [3]J_3 \,. (25) Thus, we define the inelastic potential dependent on the modified stress invariants as ω=ω(¯)=ω(I^1¯,J^2¯,J^3¯).ω=ω( )=ω ( I_1 , J_2 , J_3 )\,. (26) Unlike the previously mentioned condition for the elastic potential, the inelastic potential must be convex with respect to the modified stress invariants for the dissipation inequality to be fulfilled. Noteworthy, we could enhance the existing arguments as long as the convexity condition is satisfied. For instance, to increase network flexibility and expressibility, the negative stress invariants, or additional invariants, e.g., the principal invariants [Holthusen et al., 2026], can be appended to the argument of the inelastic potential network. These alternative constructions are presented in Appendix A.3. Following the approach of the previously presented general input-convex KAN, a convex, zero-valued and nonnegative inelastic potential ω is constructed by applying ℋH-operator (Equation 20) as a postprocessing activation on the output of the monotonic input-convex KAN, ωKANω^KAN, i.e., ω=ℋ(ωKAN)=ωKAN−ωKAN|¯=−[∂ω/∂I^1¯∂ω/∂J^2¯∂ω/∂J^3¯]|¯=0⋅[I^1¯J^2¯J^3¯].ω=H(ω^KAN)=ω^KAN- .ω^KAN |_ =0- . bmatrix∂ω/∂ I_1 \\ ∂ω/∂ J_2 \\ ∂ω/∂ J_3 bmatrix |_ =0· bmatrix I_1 \\ J_2 \\ J_3 bmatrix\,. (27) As previously stated, a further postprocessing step shown in Equation 21 can be applied to ω. Due to the special construction of the inelastic potential, which allows for a zero interval of the derivative of inelastic potential, meaning that no inelastic response is activated until a certain stress level is reached, and one single iCKAN model is sufficient to represent a viscoelastic material under relaxation. This can be represented as the rheological model in Figure 4. The elastic potential ψ is associated with the spring, and the inelastic potential ω with the dashpot. However, at this end, the training of this parameter c is unstable, leading us to choose to leave out this particular step and combine two iCKANs in parallel to express a viscoelastic material. ψ(¯e)ψ( C_e)ω(¯)ω( ) Figure 4: Representation of the iCKAN formulation in the form of a Maxwell model. The elastic potential ψ is associated with the spring, and the inelastic potential ω with the dashpot. 3.2.3 Feature dependency The appropriate functional form of the elastic and inelastic potential generally depend on material descriptors such as composition, microstructure, or processing conditions. We collect this information in a feature vector f and augment the arguments of both the elastic and inelastic potentials as ψ=ψ(I^1¯e,I^2¯e,I^3¯e,)andω=ω(I^1¯,J^2¯,J^3¯,).ψ=ψ( I_1 C_e, I_2 C_e, I_3 C_e,f) ω=ω ( I_1 , J_2 , J_3 ,f )\,. (28) This feature is solely non-mechanical and neither the elastic potential nor the inelastic potential is convex with respect to it. 3.2.4 Numerical implementation The evolution equation defined by the inelastic rate tensor, i=∂ω/∂D_i=∂ω/∂ (see also (7)), is solved numerically. The exponential integrator map has been shown to be particularly effective for inelastic problems and in the finite strain regime [Holthusen et al., 2026, 2025]. Numerically, this evolution equation can be integrated using either explicit or implicit time integration schemes. In the explicit formulation, the current inelastic stretch tensor i,tU_i,t depends solely on the previous time step and the current increment. Conversely, the implicit formulation introduces a nonlinear dependence on both current and previous states, necessitating the solution of a nonlinear system at each time step. Although the implicit scheme is generally less sensitive to time step size, it is computationally more demanding than the explicit approach. To accelerate the solution of the nonlinear system in the implicit scheme, a Liquid Time Constant Network can be employed as an alternative to the traditional iterative Newton solver [Holthusen and Kuhl, 2026]. Section 5 presents a side-by-side comparison of both formulations using synthetic data. Given its superior robustness and performance across all numerical examples discussed in the following section, this work primarily focuses on explicit time integration. Detailed information on the implicit time integration scheme is provided in Appendix A.6. Explicit time integration. Following the explicit integration scheme, the inelastic right Cauchy-Green tensor at the current timestep, i,tC_i,t, is approximated using the inelastic strain rate from the previous timestep, i,t−1D_i,t-1, as i,t=i,t−1exp(2Δti,t−1)i,t−1,i,t=i,t1/2,C_i,t=U_i,t-1 (2\, t\,D_i,t-1)\,U_i,t-1, _i,t=C_i,t^1/2\,, (29) where Δt t is the time increment between timestep t−1t-1 and t [Holthusen et al., 2026]. The overall architecture of the explicit iCKAN is illustrated in Figure 5. Noteworthy, while a Tyler expansion might be sufficient to generate the matrix squre-root under moderate deformations, a closed-form representation should be utilized for nonlinear finite strains [Holthusen et al., 2026; Hudobivnik and Korelc, 2016]. The explicit scheme proves to be computationally efficient, as no iterative solver is required. Stability is controlled through the selection of the time increment Δt t, offering a straightforward and transparent mechanism for maintaining numerical robustness. sections/standalone/iCKAN_explicit Figure 5: Explicit iCKAN architecture at timestep t. The inputs at each step are the current deformation gradient and time increment (t,Δt)(F_t, t) together with the state variables (t−1,i,t−1)(C_t-1,U_i,t-1) from the previous step. The KAN models for the elastic potential ψ and the inelastic potential ω are evaluated to update the state variables according to Equation 29. The updated state is then propagated to the next time step. In addition, the output stress P is computed by evaluating the elastic potential ψ with the updated state variables. Initial grid range. Since the KAN activation functions are represented by B-splines defined on fixed domains, suitable input grid range must be specified for each input dimension. In the present recurrent formulation, these effective inputs are not known a priori, as they are generated during the forward pass and evolve throughout training. We therefore adopt a data-driven initialization strategy to estimate appropriate grid ranges from physically meaningful tensor invariants computed on the training data. The detailed construction and its justification are provided in Appendix A.4. Further details of computation refinements, i.e. the modified definitions of square root and cubic root and the numerical implementation of the matrix square roots, are summarized in the Appendix A.5. 4 Symbolic constitutive modeling To obtain an interpretable closed-form expression of the constitutive models, the trained KANs for the elastic and inelastic potentials can be symbolified after model training. In this section, we briefly review the sparsification and symbolification process, for more details, please refer to [Abdolazizi et al., 2025; Liu et al., 2024]. Sparsification and pruning. To achieve a sparse network structure, L1 regularization is applied to the trainable parameters, which in this case correspond to the control points of the B-Spline activations. Furthermore, pruning methods can be applied to the trained KANs to remove less important connections or nodes, further simplifying the model. This can help improve interpretability by focusing on the most relevant features and reducing the overall complexity of the network. Symbolification of convex, monotonic activations. The trained KANs can be symbolified into closed-form mathematical expressions using symbolic regression techniques. The symbolification process aims to find a mathematical expression that closely approximates the behavior of the trained KANs while being more interpretable. This is achieved by searching through a space of mathematical expressions fl,j,if_l,j,i and selecting the one that best fits the data generated by the KANs. For each activation ϕl,j,i _l,j,i, we aim to find the optimal symbolically approximated activation, yl,j,i=cl,j,i⋅(fl,j,i(al,j,i⋅xl,i+bl,j,i))+dl,j,i,y_l,j,i=c_l,j,i·(f_l,j,i(a_l,j,i· x_l,i+b_l,j,i))+d_l,j,i\,, (30) so that the squared error between ϕl,j,i _l,j,i and yl,j,iy_l,j,i is minimized. For this, a library of designed symbolic functions is provided and additional affine parameters (a,b,c,d)l,j,i(a,b,c,d)_l,j,i are introduced. Furthermore, the parameters (a,c)l,j,i(a,c)_l,j,i are constrained to be nonnegative to preserve the monotonicity and convexity of the activation functions. Due to the recurrent network architecture, only the activation values from the final forward pass (i.e., the last time step) are retained. Consequently, a forward pass over the entire dataset is required to recover activation values across the full input domain. Otherwise, the extracted symbolic function does not faithfully reproduce the original B-spline activation. 5 Numerical examples In this section, we validate the proposed iCKAN framework on synthetic and experimental datasets, including VHB 4910 [Hossain et al., 2012] and temperature-dependent VHB 4905 polymer [Liao et al., 2020]. These examples demonstrate that iCKAN capture inelastic material behaviors and additional feature effects. Symbolic expressions of the trained iCKANs are presented to highlight the interpretability of the approach. For all training processes, the initial weights and biases of the KANs are initialized using a uniform random distribution. The loss function is defined as the normalized mean squared error (NMSE) between the predicted and target first Piola-Kirchhoff stresses, P and P, respectively. The error is normalized by the largest absolute entry of the target stress for each batch sample, computed over all time steps and stress components. Formally, the loss is given by ℒstress=1BTS∑b=1B∑t=1T∑s=1S|Pb,t,s−P^b,t,smaxt,s|P^b,t,s||2,L_stress= 1BTS _b=1^B _t=1^T _s=1^S | P_b,t,s- P_b,t,s _t,s | P_b,t,s | |^2\,, (31) where B is the batch size, T is the number of time steps, and S is the number of stress components. The AMSGRAD optimizer [Reddi et al., 2019], which combines the benefits of ADAM and RMSPROP, is utilized for training in combination with a cyclic learning rate scheduler [Smith, 2017], where the learning rate cyclically varies between a base learning rate and a maximum learning rate over a given step size. Furthermore, L1 regularization is applied to achieve a sparse network and gradient clipping is applied to prevent gradient explosion. 5.1 Verification on synthetic data Synthetic data are generated using an explicit time integration scheme on an iCKAN, where the corresponding KANs are replaced by analytical formulas. A Neo-Hookean free energy ψ=0.5⋅I^1¯e+0.5⋅I^3¯eψ=0.5· I_1 C_e+0.5· I_3 C_e is assumed for the elastic potential, and the inelastic potential is defined as ω=0.1⋅(I^1¯)2+(J^2¯)2ω=0.1·( I_1 )^2+( J_2 )^2. We only change the first entry of the deformation gradient from the undeformed state and apply a cyclic loading case of tension, relaxation, compression and relaxation. We set the maximum stretches to F11max=1.1,1.2,1.3F_11 =\1.1,1.2,1.3\ and set the time of tensile loading to tload=0.2,0.5,1.0t_load=\0.2,0.5,1.0\ s, resulting into different strain rates. At the end, the dataset consists of deformation gradients F and time increments Δt t as inputs and the corresponding first Piola–Kirchhoff stresses P as outputs. Model training. For the training, we use only the tensile loading and relaxation part of the datasets with tload=0.2,0.5t_load=\0.2,0.5\ s for all three stretch levels. We consider an iCKAN which can be illustrated as the rheological model in Figure 4. Both the elastic and inelastic potentials are modeled by input-convex KANs with topology [3,1]. The initial grid ranges are approximated as described in Appendix A.4. To improve robustness and account for previously unseen input ranges, the grid bounds are extended by 20%. For arguments that may attain negative values, the lower bound is chosen symmetrically with respect to the upper bound. We examine both the explicit and implicit variant of iCKAN. Training is performed with an initial learning rate of 5⋅10−45· 10^-4 and a cyclic scheduler between 5⋅10−55· 10^-5 and 1⋅10−41· 10^-4 with a step size of 50 iterations. L1 regularization of 10−510^-5 and gradient clipping at 0.1 are applied. For the implicit iCKAN, the evolution penalty with λevo=1000 _evo=1000 and increased by a factor of two every 250 epochs. Symbolification. We demonstrate the symbolification on the trained iCKAN. All activations in the KAN for elastic potential are approximated using a library of convex, non-decreasing functions f(x)f(x) on [0,∞)[0,∞). This process is illustrated in LABEL:subfig-1:psi and the full set of candidate functions is listed in Appendix A.7. The symbolic formulation of the output inelastic potential ωKANω^KAN can be obtained analogously while the enforcement of convexity, zero-valued at its origin and non-negativity is achieved by applying ℋH-operator in Equation 27 on the symbolified ωKANω^KAN. However, for the presented single-layer KAN, the network output reduces to ωKAN=∑i=1n0ϕ0,j,i(x0,i),ω^KAN= _i=1^n_0 _0,j,i(x_0,i)\,, (32) In this case, it is equivalent to apply ℋH-operator to the full output ωKANω^KAN or to each activation function ϕ0,0,i _0,0,i, ℋ(ωKAN)=∑i=1n0ℋ(ϕ0,0,i(x0,i)).H(ω^KAN)= _i=1^n_0H( _0,0,i(x_0,i))\,. (33) Consequently, the desired inelastic potential can be obtained by directly symbolizing ℋ(ϕ0,j0,i)H( _0,j0,i) using candidate functions that are convex, nonnegative, and zero at the origin (see Appendix A.7). The resulting symbolic expression inherently satisfies all constraints without requiring further postprocessing, i.e., ω=ωKANω=ω^KAN. This procedure is illustrated in LABEL:subfig-2:omega. However, this simplification does not apply to deeper network architectures. The final symbolic expressions for the elastic and inelastic potential functions identified by the implicit and explicit iCKAN approaches are presented in Figure 7c. The corresponding predictions are shown in Figure 7a and Figure 7b. More detailed results regarding model validation with synthetic data are provided in Appendix A.8. The predictions generated by iCKAN align closely with both the training and test datasets. Interestingly, both time integration approaches produced comparable results on this simple synthetic dataset, with only minor differences in the discovered symbolic forms of the elastic and inelastic potentials. This demonstrates that, in principle, either approach leads to similar outcomes and can be effectively employed in practice to discover material models using iCKANs. overpic[width=54.62706pt]figures/KAN_31.png (-1.0,28.0) [width=16.84518pt]figures/spline/psi/sp_0_0_0.png (29.0,28.0) [width=16.84518pt]figures/spline/psi/sp_0_1_0.png (58.0,28.0) [width=16.84518pt]figures/spline/psi/sp_0_2_0.png (2.0,-18.0)$ I_1 C_e$ (32.0,-18.0)$ I_2 C_e$ (63.0,-18.0)$ I_3 C_e$ (32.0,107.0)$ψ^KAN$ overpic ⟶symb. symb. overpic[width=54.62706pt]figures/KAN_31.png (0.0,28.0) [width=16.84518pt]figures/spline_symb/psi/sp_0_0_0.png (29.0,28.0) [width=16.84518pt]figures/spline_symb/psi/sp_0_1_0.png (58.0,28.0) [width=16.84518pt]figures/spline_symb/psi/sp_0_2_0.png (2.0,-18.0)$ I_1 C_e$ (32.0,-18.0)$ I_2 C_e$ (63.0,-18.0)$ I_3 C_e$ (32.0,107.0)$ψ^KAN$ overpic (a) overpic[width=54.62706pt]figures/KAN_31.png (0.0,28.0) [width=16.84518pt]figures/spline/omega/sp_0_0_0.png (29.0,28.0) [width=16.84518pt]figures/spline/omega/sp_0_1_0.png (58.0,28.0) [width=16.84518pt]figures/spline/omega/sp_0_2_0.png (2.0,-18.0)$ I_1 $ (28.0,-18.0)$ J_2 $ (60.0,-18.0)$ J_3 $ (32.0,107.0)$ω^KAN$ overpic ⟶symb. symb. overpic[width=54.62706pt]figures/KAN_31.png (0.0,28.0) [width=16.84518pt]figures/spline_symb/omega/sp_0_0_0.png (29.0,28.0) [width=16.84518pt]figures/spline_symb/omega/sp_0_1_0.png (59.0,28.0) [width=16.84518pt]figures/spline_symb/omega/sp_0_2_0.png (2.0,-18.0)$ I_1 $ (28.0,-18.0)$ J_2 $ (60.0,-18.0)$ J_3 $ (32.0,107.0)$ω^KAN$ overpic (b) Figure 6: Symbolification of the activation functions of the trained iCKAN on synthetic data. Black curves denote the learned B-spline activations, and red curves their symbolified approximations. (a) KAN-based elastic potential: convex, non-decreasing activations are approximated by symbolic functions that preserve convexity and monotonicity over the admissible domain. (b) KAN-based inelastic potential: activations are transformed by ℋH-operator (Equation 20) and then symbolified by convex, nonnegative functions that are zero-valued at the origin. (a) Explicit time integration(b) Implicit time integrationtt[s] Elastic potential ψ(¯e)=a⋅I^1¯e+b⋅I^2¯e+c⋅I^3¯eψ( C_e)=a· I_1 C_e+b· I_2 C_e+c· I_3 C_e Material parameters Explicit Implicit a b c a b c 0.0 0.2194 0.5027 0.2722 0.1093 0.5255 Inelastic potential ω(¯)=a⋅(I^1¯)2+b⋅(J^2¯)2+c⋅(J^3¯)2ω( )=a·( I_1 )^2+b·( J_2 )^2+c·( J_3 )^2 Material parameters Explicit Implicit a b c a b c 0.109 0.9394 0.4286 0.1045 1.2613 0.6397 (c) Discovered elastic and inelastic potential functions for the synthetic dataset Figure 7: Results of the symbolified iCKAN on synthetic compressible material data using (a) the explicit time integration scheme and (b) the implicit time integration scheme, using the expressions presented in Table (c). The gray regions in (a) and (b) indicate the range of the training data, and the colored dots represent the corresponding reference data. 5.2 Model discovery for experimental data To access the practical applicability of iCKANs, we examine the framework on experimental data of viscoelastic polymers. Due to the fact that the provided experimental data are one-dimensional, the information of volumetric modes are lacking. Thus, we assume incompressible materials [Holthusen et al., 2026; Abdolazizi et al., 2024]. The detailed hyperparameters are provided in Appendix A.9, together with the loss convergence curves over the training epochs. 5.2.1 Viscoelastic model discovery of VHB 4910 polymer Very-High-Bond (VHB) 4910 is a polymer exhibiting highly nonlinear and viscoelastic behavior. Experimental uniaxial loading–unloading data for VHB 4910 are reported by Hossain et al. [2012] for four maximum stretch levels, F11=1.5,,2.0,,2.5,,3.0F_11=\1.5,,2.0,,2.5,,3.0\, tested under three different stretch rates, F˙11=0.01,,0.03,,0.05s−1 F_11=\0.01,,0.03,,0.05\\,s^-1. In our framework, rate dependence is incorporated directly through the network input pair (,Δt)(F, t), allowing the model to learn time-dependent behavior without introducing explicit viscosity terms. To capture more complex material behaviors, we utilize two iCKANs in parallel, motivated by classical rheological representations of viscoelastic materials, as illustrated in Figure 8b. In continuum mechanics, the total stress response of viscoelastic solids is often represented as the superposition of multiple parallel mechanisms, each corresponding to a distinct elastic potential and inelastic potential. The two branches of the parallel structure are indexed by superscripts (∙)1^1( ) and (∙)2^2( ), which share the same input deformation gradient F. Each branch learns its own elastic potential ψi(i¯e)^iψ(^i C_e) and inelastic potential ωi(i¯)^iω(^i ), while enforcing thermodynamic consistency. The stress contribution i^iP of each branch is obtained consistently via differentiation of the learned potentials. The total stress response is then computed as the additive superposition =1+2P=^1P+^2P, which is consistent with classical rheological network models and enables the representation of multiple concurrent elastic and inelastic mechanisms with different characteristic responses. Model training. For training, we use the dataset corresponding to the highest stretch level, F11=3.0F_11=3.0, while the remaining stretch levels are reserved for validation to assess generalization for unseen data. We approximate the initial grid range as previously and extend them by 20% on both sides. overpic[height=455.24408pt]results/VHB4910/KAN_321_psi1.png (-1.0,15.0) [height=68.28383pt]results/VHB4910/spline_psi1/sp_0_0_0.png (30.0,15.0) [height=68.28383pt]results/VHB4910/spline_psi1/sp_1_0_0.png (15.0,62.0) [height=68.28383pt]results/VHB4910/spline_psi1/sp_0_1_0.png (10.0,-14.0)$ I_1^^1 C_e$ (42.0,-14.0)$ I_2^^1 C_e$ (74.0,-14.0)$ I_3^^1 C_e$ (38.0,107.0)$^1ψ^KAN$ overpic (a) ψ1(1¯e)^1ψ(^1 C_e) overpic[height=455.24408pt]results/VHB4910/KAN_321_omega1.png (45.0,15.0) [height=68.28383pt]results/VHB4910/spline_omega1/sp_0_2_1.png (61.0,62.0) [height=68.28383pt]results/VHB4910/spline_omega1/sp_1_1_0.png (10.0,-14.0)$ I_1^^1 $ (42.0,-14.0)$ J_2^^1 $ (74.0,-14.0)$ J_3^^1 $ (38.0,107.0)$^1ω^KAN$ overpic (b) ω1(1¯)^1ω(^1 ) overpic[height=455.24408pt]results/VHB4910/KAN_321_psi2.png (-1.0,15.0) [height=68.28383pt]results/VHB4910/spline_psi2/sp_0_0_0.png (15.0,15.0) [height=68.28383pt]results/VHB4910/spline_psi2/sp_0_0_1.png (30.0,15.0) [height=68.28383pt]results/VHB4910/spline_psi2/sp_0_1_0.png (45.0,15.0) [height=68.28383pt]results/VHB4910/spline_psi2/sp_0_1_1.png (15.0,62.0) [height=68.28383pt]results/VHB4910/spline_psi2/sp_1_0_0.png (61.0,62.0) [height=68.28383pt]results/VHB4910/spline_psi2/sp_1_1_0.png (10.0,-14.0)$ I_1^^2 C_e$ (42.0,-14.0)$ I_2^^2 C_e$ (74.0,-14.0)$ I_3^^2 C_e$ (38.0,107.0)$^2ψ^KAN$ overpic (c) ψ2(2¯e)^2ψ(^2 C_e) ψ1(1¯e)^1ψ(^1 C_e)ψ2(2¯e)^2ψ(^2 C_e)ω1(1¯)^1ω(^1 )ω2(2¯)^2ω(^2 )=P=1+2^1P+^2P (a) Discovered model for VHB4910 polymer(b) Rheological model Figure 8: KAN architectures of discovered model for VHB4910 polymer. (b) Rheological model consisting two iCKANs in parallel combination. Both iCKAN branches receive the same input deformation gradient F and each learns an elastic potential and an inelastic potential. The first branch identifies ψ1(1¯e)^1ψ(^1 C_e) and ω1(1¯)^1ω(^1 ) and produces the stress 1^1P while preserving thermodynamic consistency. The second branch identifies ψ2(2¯e)^2ψ(^2 C_e) and ω2(2¯)^2ω(^2 ) and produces the stress 2^2P under the same thermodynamic constraints. The total stress response is obtained by summing the stresses of both branches, =1+2P=^1P+^2P. Symbolification and results. At the end, we utilize the same library of convex, monotoonic functions as in the previous section to symbolify the activations of the trained model. The symbolified discovered KAN models of the corresponding elastic and inelastic potentials are shown in Figure 8a with their corresponding symbolic functions in Figure 9b. The predicted results of the symbolified iCKAN model are shown in Figure 9a. TrainingTestingTestingTesting(a) Predictions of the symbolified iCKANs for VHB 4910 Elastic potential ψ1(1¯e)=a⋅exp(b⋅I^1¯e1)+c⋅I^2¯e1^1ψ(^1 C_e)=a· (b· I_1^^1 C_e)+c· I_2^^1 C_e Material parameters a b c 34.654 0.1204 3.7902 Inelastic potential ω1(1¯)=ℋ(a⋅exp(b⋅exp(c⋅J^2¯1)))^1ω(^1 )=H(a· (b· (c· J_2^^1 ))) Material parameters a b c 0.0195 0.0002 0.036 Elastic potential ψ2(2¯e)=a⋅I^1¯e2+b⋅I^2¯e2+c⋅exp(d⋅I^1¯e2)^2ψ(^2 C_e)=a· I_1^^2 C_e+b· I_2^^2 C_e+c· (d· I_1^^2 C_e) Material parameters a b c d 2.5911 1.8997 68.73 0.0299 (b) Discovered elastic and inelastic potential functions for VHB 4910 polymer Figure 9: Discovered model for VHB4910 polymer. (a) Predictions by the symbolified trained iCKAN under three different constant loading/unloading rates, denoted as F˙11 F_11. The prediction of the trained iCKAN model is represented by the solid lines and the experimentally measured data is represented by dots in corresponding colors. Only the experimental data corresponding to F11max=3.0F_11^max=3.0 [-] are utilized for training. (b) listes the corresponding symbolic functions. 5.2.2 Thermo-viscoelastic model discovery of VHB 4905 polymer Building on the previous section, we next consider the identification of material models augmented by an additional feature vector f. Our goal is to find a single iCKAN expression capable of describing the material behavior under varying non-mechanical conditions. To this end, we study the VHB 4905 polymer, which exhibits highly deformable, viscoelastic, and temperature-sensitive behavior. Uniaxial experimental data for VHB 4905 polymer at different temperatures, θ=0,10,20,40,60,80θ=\0,10,20,40,60,80\ ∘C, under multiple strain rates, F11=0.03,0.05,0.1F_11=\0.03,0.05,0.1\ s-1, and maximum stretch ratios, F11=2.0,3.5,4.0F_11=\2.0,3.5,4.0\, are reported in [Liao et al., 2020]. Temperature as feature. The strain rate and deformation information are naturally incorporated into the model inputs. In addition, the temperature θ is explicitly included as a feature in the arguments of the elastic potential, i.e., [Abdolazizi et al., 2025] ψi(i¯e)=ψi(I^1¯ei,I^2¯ei,I^3¯ei,θ).^iψ(^i C_e)=^iψ( I_1^^i C_e, I_2^^i C_e, I_3^^i C_e,θ)\,. (34) It is important to note that this approach differs from classical thermo-mechanical constitutive modeling, where temperature is included via either thermal strain decomposition or temperature-dependent material parameters. Instead, this feature-augmented construction provides a general, data-driven mechanism to capture the influence of unknown non-mechanical factors on the elastic potential, and consequently, on the material stress response. The initial grid range for this temperature feature is set directly by the experimental temperature range, [0,80][0,80] ∘C, while the grid ranges for the invariants are approximated as before and extend by 20% on both ends. Although it is reasonable to assume the feature has influence on both elastic and inelastic potential, we choose to omit it for the inelastic potential due to network robustness. Symbolification. We symbolify the trained activations, except for those corresponding to the temperature feature in the first layer of each KAN for elastic energy, using the same library of convex, monotonic functions as in the previous section. The resulting KAN architectures are shown in Figure 10. For each branch, there remains an activation function associated with temperature in the elastic potential ψi(i¯e)^iψ(^i C_e), which we denote as g1(θ)^1g(θ) and g2(θ)^2g(θ) for the first and second branch, respectively. overpic[height=91.04742pt]results/VHB4905/KAN_421_psi1.png (-1.0,12.0) [height=13.65675pt]results/VHB4905/spline_psi1/sp_0_0_0.png (25.0,12.0) [height=13.65675pt]results/VHB4905/spline_psi1/sp_0_1_0.png (77.0,12.0) [height=13.65675pt]results/VHB4905/spline_psi1/sp_0_3_0.png (25.0,52.0) [height=13.65675pt]results/VHB4905/spline_psi1/sp_1_0_0.png (10.0,-13.0)$ I_1^^1 C_e$ (33.0,-13.0)$ I_2^^1 C_e$ (60.0,-13.0)$ I_3^^1 C_e$ (88.0,-13.0)$θ$ (45.0,88.0)$^1ψ^KAN$ overpic (a) overpic[height=91.04742pt]results/VHB4905/KAN_321_omega1_VHB4905.png (30.0,15.0) [height=13.65675pt]results/VHB4905/spline_omega1/sp_0_1_0.png (62.0,15.0) [height=13.65675pt]results/VHB4905/spline_omega1/sp_0_2_0.png (14.5,62.0) [height=13.65675pt]results/VHB4905/spline_omega1/sp_1_0_0.png (10.0,-15.0)$ I_1^^1 $ (42.0,-15.0)$ J_2^^1 $ (74.0,-15.0)$ J_3^^1 $ (38.0,105.0)$^1ω^KAN$ overpic (b) overpic[height=91.04742pt]results/VHB4905/KAN_421_psi2.png (38.0,12.0) [height=13.65675pt]results/VHB4905/spline_psi2/sp_0_1_1.png (65.0,12.0) [height=13.65675pt]results/VHB4905/spline_psi2/sp_0_2_1.png (90.0,12.0) [height=13.65675pt]results/VHB4905/spline_psi2/sp_0_3_1.png (65.0,52.0) [height=13.65675pt]results/VHB4905/spline_psi2/sp_1_1_0.png (10.0,-13.0)$ I_1^^2 C_e$ (33.0,-13.0)$ I_2^^2 C_e$ (60.0,-13.0)$ I_3^^2 C_e$ (88.0,-13.0)$θ$ (45.0,88.0)$^2ψ^KAN$ overpic (c) overpic[height=91.04742pt]results/VHB4905/KAN_321_omega2_VHB4905.png (-1.0,15.0) [height=13.65675pt]results/VHB4905/spline_omega2/sp_0_0_0.png (30.0,15.0) [height=13.65675pt]results/VHB4905/spline_omega2/sp_0_1_0.png (62.0,15.0) [height=13.65675pt]results/VHB4905/spline_omega2/sp_0_2_0.png (14.5,62.0) [height=13.65675pt]results/VHB4905/spline_omega2/sp_1_0_0.png (10.0,-15.0)$ I_1^^2 $ (42.0,-15.0)$ J_2^^2 $ (74.0,-15.0)$ J_3^^2 $ (38.0,105.0)$^2ω^KAN$ overpic (d) Figure 10: Discovered model for VHB4905 thermo polymer As these activations are not constrained to be monotonic or convex, their shapes can be arbitrary. To capture this flexibility, we adopt higher-order polynomial functions. This approach allows us to automatically discover a symbolic expression for g1(θ)^1g(θ). For g2(θ)^2g(θ), however, a single symbolifc function did not provide sastisfactory fit. Therfore, we allow for piecewise symbolification, and manually approximate it using two piecewise quadratic functions defined over the temperature intervals [0,40)[0,40) and [40,80][40,80] ∘C, respectively. The coefficients are chosen such that C1C^1-continuity is enforced at θ=40θ=40 ∘C. Results. The predicted results of the symbolified iCKAN model are shown in Figure 11. The symbolified formula for the discoverd model are listed in Figure 12. (a) Training data (b) Testing data Figure 11: (a) Training set and (b) validation set of discovered model for the experimental data of VHB 4905 polymer under three different constant loading/unloading rates, denoted as F˙11 F_11. The prediction of the trained iCKAN model is represented by the solid lines and the experimentally measured data is represented by the dotts in the corresponding colors. The iCKAN framework successfully captures the temperature-dependent viscoelastic behavior of VHB 4905. Symbolification of the learned activations provides interpretable functional forms, revealing how temperature influences the elastic potentials in each branch. Overall, this demonstrates the model’s ability to combine predictive accuracy with physical interpretability under varying environmental conditions. Elastic potential ψ1(1¯e)=a⋅exp(b⋅I^1¯e1+c⋅I^2¯e1+¹g(θ))^1ψ(^1 C_e)=a· (b· I_1^^1 C_e+c· I_2^^1 C_e+¹g(θ)) Material parameters a b c 0.0349 0.0855 0.0864 Temperature dependency g1(θ)=10.0746⋅(0.5239−0.0111⋅θ)3^1g(θ)=10.0746·(0.5239-0.0111·θ)^3 Inelastic potential ω1(1¯)=ℋ(a⋅exp(b⋅exp(c⋅J^2¯1)+d⋅exp(e⋅J^3¯1)))^1ω(^1 )=H(a· (b· (c· J_2^^1 )+d· (e· J_3^^1 ))) Material parameters a b c d e 31.918 0.1077 4.7 0.0382 10.0 Elastic potential ψ2(2¯e)=a⋅exp(b⋅I^2¯e2+c⋅exp(d⋅I^3¯e2)+2g(θ))^2ψ(^2 C_e)=a· (b· I_2^^2 C_e+c· (d· I_3^^2 C_e)+^2g(θ)) Material parameters a b c d 0.319 0.0261 0.0921 10.0 Temperature dependency g2(θ)=−1.7498⋅10−4⋅θ2+3.256⋅10−3⋅θ−1.8712⋅10−1,0≤θ<40−3.51⋅10−5⋅θ2+2.344⋅10−3⋅θ−7.9019⋅10−2,40≤θ≤80^2g(θ)= cases-1.7498· 10^-4·θ^2+3.256· 10^-3·θ-1.8712· 10^-1,&0≤θ<40\\[6.0pt] -3.51· 10^-5·θ^2+2.344· 10^-3·θ-7.9019· 10^-2,&40≤θ≤ 80 cases Inelastic potential ω2(2¯)=ℋ(a⋅exp(b⋅I^1¯2+c⋅exp(d⋅J^2¯2)))^2ω(^2 )=H(a· (b· I_1^^2 +c· (d· J_2^^2 ))) Material parameters a b c d 0.0589 0.2543 0.0264 2.8013 Figure 12: Discovered elastic and inelastic potential functions for VHB 4905 thermo polymer. 6 Discussion Developing accurate and interpretable constitutive models for materials undergoing nonlinear inelastic deformation, particularly at finite strains, remains a central challenge in computational mechanics. Although recent advances in data-driven modeling have accelerated model discovery, a persistent trade-off exists between predictive accuracy and physical interpretability. Thermodynamically consistent approaches often rely on highly flexible neural networks, which can lack transparency, or on symbolic methods with limited expressiveness [Holthusen et al., 2026; Flaschel et al., 2023]. This work bridges this divide by introducing inelastic Constitutive Kolmogorov–Arnold Networks (iCKANs), which combine a generalized inelastic constitutive framework with the interpretability of Kolmogorov–Arnold Networks (KANs). This integration enables iCKANs to achieve predictive accuracy, interpretability, and robust extrapolation simultaneously. Their flexible architecture also facilitates symbolic extraction of learned potentials. Building upon previous KAN-based frameworks for hyperelastic materials, iCKANs establish a new paradigm for interpretable, data-driven discovery of inelastic material models [Abdolazizi et al., 2025; Thakolkaran et al., 2025]. Interpretability. The iCKAN framework shares the constitutive structure of inelastic Constitutive Neural Networks (iCANNs) [Holthusen et al., 2026], but differs fundamentally in its representation of the elastic and inelastic potentials. Whereas iCANNs employ custom-built convex neural networks, iCKANs use input-convex KANs, enabling a more compact and expressive formulation. The use of trainable B-spline activations enables iCKANs to achieve high expressiveness with substantially reduced network complexity. A major advantage of KANs is their built-in support for activation-level symbolic regression, which facilitates the direct extraction of interpretable, closed-form analytical expressions for both elastic and inelastic potentials. This capability preserves flexibility in data-driven discovery while yielding transparent, physically meaningful models with strong predictive performance [Abdolazizi et al., 2025]. Furthermore, iCKANs naturally accommodate additional non-mechanical features, such as temperature, whose effects can be extracted as explicit mathematical expressions, providing valuable physical insight. Extrapolability. The combination of a thermodynamically consistent constitutive framework and convex-network architecture ensures that the learned elastic and inelastic potentials satisfy essential physical constraints across their entire input domain. Importantly, these constraints are retained in the closed-form symbolic expressions obtained through the KAN symbolification process. This structure enhances robustness and enables physically consistent extrapolation beyond the training regime [Abdolazizi et al., 2025]. As demonstrated in synthetic benchmarks (Section 5.1), models trained on limited deformation data can still provide reliable predictions for unseen, larger deformations using both the spline-based KAN and its symbolic form. Computational efficiency. By promoting sparsity in the KAN architecture, the symbolic regression process yields concise, closed-form expressions for both elastic and inelastic potentials. These analytical formulations can be seamlessly integrated into finite element solvers, eliminating the need for network evaluation during simulation. As a result, the iCKAN approach introduces no additional computational overhead at deployment [Abdolazizi et al., 2025]. Choice of input representation. This study employs the principal invariants of the elastic right Cauchy–Green tensor as inputs for the elastic potential, but alternative objective representations—such as principal stretches—are equally valid in the iCKAN formulation [Abdolazizi et al., 2025]. For the inelastic potential, invariant-based stress measures are used, with the option to augment input sets by combining stress and deformation invariants. The selection of input representation influences both predictive accuracy and interpretability, and a systematic comparison across different material classes is a promising direction for future research [Holthusen et al., 2026]. Predictive performance. The results demonstrate that iCKANs deliver strong predictive performance on both synthetic and experimental viscoelastic datasets, including cases with feature-dependent material behavior. However, in certain scenarios, conventional neural network models may still yield lower prediction errors. The iCKAN framework deliberately prioritizes interpretability and physical structure, which may result in modest accuracy trade-offs. Closing this gap remains an open challenge, motivating further work on training strategies, architectural optimization, and advanced regularization techniques to enhance predictive performance while retaining interpretability. Limitations and future work. While this work primarily addresses viscoelastic material behavior, the iCKAN framework is broadly applicable to other classes of inelasticity, including plasticity, damage, growth and remodeling [Boes et al., 2026; Holthusen et al., 2026, 2025]. Extensions to anisotropic materials and structural-scale implementations are also feasible [Holthusen and Kuhl, 2026; Kalina et al., 2025; Thakolkaran et al., 2025]. These directions open new possibilities for feature-aware material characterization and data-driven design. From a numerical standpoint, KANs require the specification of spline grid ranges, which makes initialization more involved than in standard multilayer perceptrons. Although the present approximation strategy proved effective, future research should focus on developing robust adaptive grids and training techniques. Additional opportunities include investigating partially input-convex architectures, advanced pruning and convexification strategies, and automated hyperparameter selection, potentially leveraging large language models or other AI-driven tools [Hospedales et al., 2021; Tacke et al., 2025]. Efficient, standardized integration into finite element solvers for both forward and inverse boundary value problems represents another key step for broad adoption. Conclusion. Reliable and interpretable constitutive models for inelastic materials under finite deformations are essential for fields where complex, feature-dependent responses arise, such as biomedical engineering. The iCKAN framework presented in this work provides a robust, physically consistent, and data-driven solution for automated discovery of inelastic material models. By incorporating input-convex Kolmogorov–Arnold Networks, the approach achieves a unique combination of flexibility, predictive accuracy, and transparent symbolic representations, advancing the state of the art in interpretable inelastic material modeling. Appendix A Appendix A.1 Generation of a convex function from a convex, non-decreasing function Let f()f(x) be a multivariate funciton of x, which is convex and non-decreasing with respect to each component of x. We prove that f~()=f()−f(α)−∇f()|=α⋅(−α) f(x)=f(x)-f(α)- .\ (x) |_x=α·(x-α) (A.1) is convex with respect to each component of x, nonnegative and has zero value and zero-valued derivative at =αx=α. Convexity. Suppose f()f(x) is a twice-diffentiable function, it holds that the second partial derivate with respect to each argument must be nonnegative, i.e., ∂2f/∂xi2≥0∂^2f/∂ x_i^2≥ 0 for all i. The first partial derivative of f~() f(x) with respect to xix_i is given as ∂f~()∂xi=∂f()∂xi−∂f()∂xi|=α, ∂ f(x)∂ x_i= ∂ f(x)∂ x_i- . ∂ f(x)∂ x_i |_x=α\,, (A.2) and the second partial derivative is ∂2f~()∂xi2=∂2f()∂xi2≥0for all i. ∂^2 f(x)∂ x_i^2= ∂^2f(x)∂ x_i^2≥ 0 all i\,. (A.3) Therefore, f~() f(x) is convex with respect to each component of x. Nonnegativity. Since f()f(x) is convex and non-decreasing with respect to each component of x, it must hold that for any x and α, f()≥f(α)+∇f()|=α⋅(−α).f(x)≥ f(α)+ .\ (x) |_x=α·(x-α)\,. (A.4) Rearranging this inequality gives f()−f(α)−∇f()|=α⋅(−α)≥0,f(x)-f(α)- .\ (x) |_x=α·(x-α)≥ 0\,, (A.5) which implies that f~()≥0 f(x)≥ 0 for all x. Zero value and zero derivative at =αx=α. Evaluating f~() f(x) at =αx=α gives f~(α)=f(α)−f(α)−∇f()|=α⋅(α−α)=0. f(α)=f(α)-f(α)- .\ (x) |_x=α·(α-α)=0\,. (A.6) The gradient of f~() f(x) at =αx=α is ∇f~()|=α=∇f()|=α−∇f()|=α=. .\ ∇ f(x) |_x=α= .\ (x) |_x=α- .\ (x) |_x=α=0\,. (A.7) Thus, f~() f(x) has a zero value and derivative at =αx=α. A.2 Stress-free undeformed state To ensure that the undeformed state is correctly represented, that is, that there should be no energy or stress at zero elastic deformation [Holzapfel, 2000], the elastic potential function must satisfy the following conditions ψ(¯e=)=0and∂ψ∂¯e|¯e==.ψ( C_e=I)=0 . ∂ψ∂ C_e |_ C_e=I=0\,. (A.8) This can be achieved by adding a correction term for stress and potential, ψσψ^σ and ψϵψ^ε, respectively, to the output elastic potential of the KAN, ψKANψ^KAN, i.e., ψ=ψKAN+ψσ(J)+ψϵψ=ψ^KAN+ψ^σ(J)+ψ^ε. These two terms are calculated dependent on the output elastic potential at the undeformed state. For more detailed explanation, please refer to Appendix D in [Abdolazizi et al., 2025]. In our case, the derivatives of the modified invariants w.r.t. ¯e C_e are zero at the undeformed state [Thakolkaran et al., 2025], i.e., ∂I^1¯e∂¯e|¯e==∂I^2¯e∂¯e|¯e==∂I^3¯e∂¯e|¯e==. . ∂ I_1 C_e∂ C_e |_ C_e=I= . ∂ I_2 C_e∂ C_e |_ C_e=I= . ∂ I_3 C_e∂ C_e |_ C_e=I=0\,. (A.9) Therefore, the stress-free undeformed state is automatically ensured. The potential-free undeformed state is achieved by subtracting the output elastic potential evaluated at in the undeformed state, ψ=ψKAN−ψKAN|¯e=.ψ=ψ^KAN- .ψ^KAN |_ C_e=I\,. (A.10) A.3 Alternative choices of argument for inelastic potential Adding negative stress invariants. For more network expressibility, the negative invariants can be added additionally to the positive invariants to enhance the input arguments of KAN for the inelastic potential, i.e., ω=ω(I^1¯,J^2¯,J^3¯,−I^1¯,−J^2¯,−J^3¯).ω=ω ( I_1 , J_2 , J_3 ,- I_1 ,- J_2 ,- J_3 )\,. (A.11) This approach increases the expressivity of the inelastic potential network by additional activations that are monotonic decreasing and convex with respect to the invariants. In this way, the convexity with respect to the invariants is preserved. To ensure that the output inelastic potential is zero-valued at its origin and non-negative, the same transformation step in Equation 27 should be carried out. Adding principal invariants. Additional stress invariants can be added to the input arguments of the inelastic potential network to enhance the expressivity, e.g., the classical invariants I2¯I_2 and I3¯I_3 . These invariants can be expressed in terms of the stress invariants as I2¯=(I1¯)26+J2¯,I3¯=(I1¯)327+23I1¯J2¯+J3¯.I_2 = (I_1 )^26+J_2 \,, I_3 = (I_1 )^327+ 23I_1 J_2 +J_3 \,. (A.12) Analogously to the modification of the stress invariants, we take the square root and cubic root of the second and third principal invariants, respectively, and achieve the following modified invariants, I^2¯=I2¯,I^3¯=I3¯3. I_2 = I_2 \,, I_3 = [3]I_3 \,. (A.13) Thus, the expression of the inelastic potential can be reformulated as ω=ω(I1¯,J^2¯,J^3¯,I^2¯(I1¯,J2¯),I^3¯(I1¯,J2¯,J3¯)).ω=ω (I_1 , J_2 , J_3 , I_2 (I_1 ,J_2 ), I_3 (I_1 ,J_2 ,J_3 ) )\,. (A.14) The guarantee of thermodynamic consistency is prove in Holthusen et al. [2026]. Analogously, the transformation step in Equation 27 should follow. A.4 Estimation for initial grid range of iCKANs The activation functions in KANs are represented by B-splines defined on fixed grids. In this work, we adopt the spline grid initialization strategy proposed by Thakolkaran et al. [2025]. Accordingly, each input dimension requires a prescribed grid range that specifies the interval over which the activation functions are explicitly represented. Outside this interval, linear extrapolation is applied, which may lead to numerical instabilities if the inputs extend significantly beyond the grid bounds. Consequently, a reliable choice of the grid range for the first KAN layer is essential for stable training and evaluation. Although one could assume a generally sufficiently large range, this would lead to a reduction of training efficiency and resolution. In contrast to prior work focusing on hyperelasticity [Abdolazizi et al., 2025; Thakolkaran et al., 2025], the effective network inputs in the present iCKAN formulation are not directly given by the dataset (,)(F,P), but instead arise within a recurrent constitutive update and depend on intermediate model predictions. Their distribution is therefore not known a priori and cannot be determined through a simple preprocessing step. To address this difficulty, we introduce an estimation procedure for suitable grid ranges based on invariant measures computed from the input and output tensors using the training data. Specifically, the initial grid range for the KAN model representing the Helmholtz free energy ψ is defined based on the invariants of the right Cauchy-Green deformation tensor =TC=F^TF, while the initial grid range for the KAN model representing the inelastic potential ω is defined based on the invariants of the product of the right Cauchy-Green deformation tensor C and the second Piola-Kirchhoff stress tensor =−1S=F^-1P. In the following, we provide a justification for the plausibility of this estimation scheme. A.4.1 Approximation for initial grid range of KAN for inelastic potential The initial grid range of the KAN for inelastic potential should cover the range of the invariants of the symmetric elastic Mandel-like stress tensor ¯ . The elastic Mandel-like stress ¯ and the second Piola-Kirchhoff stress tensor S are defined by ¯=2¯e∂ψ∂¯eand=2i−1∂ψ∂¯ei−1, =2\, C_e ∂ψ∂ C_e =2\,U_i^-1 ∂ψ∂ C_eU_i^-1\,, (A.15) leading to the relation ¯=i−1i =U_i^-1CSU_i. We thus define ^= =CS and proove that ¯=i−1^i =U_i^-1 U_i and share the same invariants. First stress invariant. I1¯=tr(¯)=tr(i−1^i)=tr(^ii−1)=tr(^)=I1Σ^.I_1 =tr( )=tr(U_i^-1 U_i)=tr( U_iU_i^-1)=tr( )=I_1 \,. (A.16) Second stress invariant. J2¯=12tr(dev(¯)2)=12tr(¯2−23I1¯¯+19(I1¯)2)=12(tr(¯2)−59(I1¯)2)=12(tr(^2)−59(I1Σ^)2)=12tr(dev(¯)2)=J2Σ^withtr(¯2)=tr(i−1^ii−1^i)=tr(i−1^^i)=tr(^2). splitJ_2 &= 12tr(dev( )^2)= 12tr ( ^2- 23I_1 + 19(I_1 )^2I )\\ &= 12(tr( ^2)- 59(I_1 )^2)= 12(tr( ^2)- 59(I_1 )^2)\\ &= 12tr(dev( )^2)=J_2 \\ &with ( ^2)=tr(U_i^-1 U_iU_i^-1 U_i)=tr(U_i^-1 U_i)=tr( ^2)\,. split (A.17) Third stress invariant. J3¯=13tr(dev(¯)3)=13tr(¯3−I1¯¯2+13(I1¯)2¯−127(I1¯)3)=13(tr(¯3)−I1¯tr(¯2)+827(I1¯)3)=13(tr(^3)−I1Σ^tr(^2)+827(I1Σ^)3)=13tr(dev(^)3)=J3Σ^withtr(¯3)=tr(i−1^ii−1^ii−1^i)=tr(i−1^^^i)=tr(^3). splitJ_3 &= 13tr(dev( )^3)= 13tr ( ^3-I_1 ^2+ 13(I_1 )^2 - 127(I_1 )^3I )\\ &= 13(tr( ^3)-I_1 tr( ^2)+ 827(I_1 )^3)= 13(tr( ^3)-I_1 tr( ^2)+ 827(I_1 )^3)\\ &= 13tr(dev( )^3)=J_3 \\ &with ( ^3)=tr(U_i^-1 U_iU_i^-1 U_iU_i^-1 U_i)=tr(U_i^-1 U_i)=tr( ^3)\,. split (A.18) In conclusion, the approximated grid range using the invariants of ^= =CS is equivalent to the grid range using the invariants of ¯ . A.4.2 Approximation for initial grid range of KAN for elastic potential The initial grid range of the KAN for elastic potential should cover the range of the invariants of the elastic right Cauchy-Green deformation tensor ¯e C_e. We approximate the grid range using the invariants of the right Cauchy-Green deformation tensor C. For one cyclic tension test, the elastic deformation is always smaller than the total deformation, i.e., ¯e≤ C_e . Hence, the invariants of ¯e C_e are always smaller than or equal to the invariants of C, i.e., I1¯e≤I1I_1 C_e≤ I_1^C, I2¯e≤I2I_2 C_e≤ I_2^C and I3¯e≤I3I_3 C_e≤ I_3^C. Therefore, the grid range using the invariants of C always covers the grid range using the invariants of ¯e C_e. It is to note, this only holds for compressible cases. For incompressible materials, we assume this approach is sufficient. A.5 Refinement of computations for numerical stability Square root and cubic root. To enable the differentiability at the origin, the square root and cubic root of the stress invariants are modified as [Holthusen et al., 2026] x=x(x+ϵ1)1/2,x3=x(x+ϵ2)2/3,with ϵ1,ϵ2>0. x= x(x+ _1)^1/2, [3]x= x(x+ _2)^2/3, _1, _2>0\,. (A.19) Matrix square root. Closed-form formulation of the matrixs by Hudobivnik and Korelc [2016] are used for computing the matrix square roots, specifically, i=iU_i= C_i [Holthusen et al., 2026]. A.6 Implicit time integration scheme Following the implicit integration scheme, the following nonlinear equation hast to be solved at each timestep t to get the current inelastic stretch tensor i,tU_i,t, :=i,t−1−i,texp(−2Δti,t)i,t=!,r:=C_i,t-1-U_i,t (-2\, t\,D_i,t)U_i,t !=0\,, (A.20) where i,tD_i,t depends nonlinearly on i,tU_i,t [Holthusen and Kuhl, 2026]. Classically, the initial guess for i,t(0)U_i,t^(0) is chosen as the previous timestep value i,t−1U_i,t-1 and the equation is solved iteratively using, e.g., the Newton-Raphson method, until the residual r is sufficiently small. In comparison to the explicit time integration scheme, this approach is more stable but however time consuming. A.6.1 Helper network for solving implicit evolution equation To overcome the drawbacks of the implicit scheme and accelerate the training, the iterative solver can be replaced by a helper neural network fN_f to directly predict i,tU_i,t [As’ ad and Farhat, 2023; Rosenkranz et al., 2024]. This helper network can be represented by a Liquid Time-Constant Network (LTC) [Hasani et al., 2021] to predict the current inelastic stretch [Holthusen and Kuhl, 2026]. For instance, we need to compute a trial value of i,tD_i,t, i,ttrialD_i,t^trial, which is iD_i evaluated with the current tC_t but the last i,t−1U_i,t-1. Note that the outcome of the network generally not statisfy the symmetric positive definite condition of i,tU_i,t. To enforce this condition, it is of advantage to work with the lower triangular matrix iL_i of the Cholesky decomposition of iU_i, i.e., i=iiTU_i=L_iL_i^T [Benoit, 1924]. The helper network is then designed to predict i,tL_i,t as i,t=f(i,t−1,t,¯i,ttrial),L_i,t=N_f(L_i,t-1,C_t, D_i,t^trial)\,, (A.21) following with i,t=i,ti,tTU_i,t=L_i,tL_i,t^T. To evaluate the accuracy of the predicted i,tU_i,t, the loss function is enhanced by a term corresponding to the evolution equation, namely the residual of the implicit scheme in Equation A.20, ℒtotal=ℒstress+λevo⋅ℒevowithℒevo=MSE()L_total=L_stress+ _evo·L_evo _evo=MSE(r) (A.22) with λevo _evo being a hyperparameter to balance the two loss parts and is heuristically set to 1e3 to 1e5 to bring the two loss parts to a similar scale. The overall architecture of the implicit iCKAN with a helper LTC is illustrated in Figure A.1. The algorithm for the implicit iCKANs with helper network is summarized in Appendix A.6. sections/standalone/iCKAN_implicit Figure A.1: Implicit iCKAN architecture with helper LTC to solve the evolution equation at timestep t. The inputs at each step are the current deformation gradient and time increment (t,Δt)(F_t, t) together with the state variables (t−1,i,t−1)(C_t-1,U_i,t-1) from the previous step. The updated state variables is achieved by solving Equation A.20 using a helper LTC network. The updated state is then propagated to the next time step. In addition, the output stress P is computed by evaluating the elastic potential function ψ with the updated state variables. Liquid time constant networks. LTCs are a type of recurrent neural networks that is designed to model temporal dependencies in sequential data. A LTC is based on the following ordinary differential equation (ODE) ddt=f(,)−α(,)(t), dhdt=f(x,h)-α(x,h)\,h(t)\,, (A.23) where h is the hidden state, x is the input, and α is a learnable parameter that controls the decay rate of the hidden state. Consequently, a LTC consists of two networks, the source network fN_f to learn the update function f, and a second network αN_α to learn the state-dependent time constant α, which is positive and scaled by the timestep Δt t. This parameter controls the dynamics of the network, specifically, the larger α, the faster the state changes. Thus, the original forward Euler scheme t+1=t+Δt⋅f(t,)h_t+1=h_t+ t· f(h_t,x) (A.24) is modified into a semi-explicit scheme, which reads t+1=t+Δt⋅(f(t,)−αtt).h_t+1=h_t+ t· (f(h_t,x)- _t\,h_t )\,. (A.25) In the present framework, the current hidden state is the previous lower Cholesky decompsition of the inelastic stretch tensor i,t−1L_i,t-1 and the current state is the current right Cauchy-Green tensor tC_t and the trial inelastic strain rate i,ttrialD_i,t^trial, i.e., i,t=i,t−1+Δt⋅(f(i,t−1,t,¯i,ttrial)−α(i,t−1,t,¯i,ttrial)i,t−1).L_i,t=L_i,t-1+ t· (N_f(L_i,t-1,C_t, D_i,t^trial)-N_α(L_i,t-1,C_t, D_i,t^trial)L_i,t-1 )\,. (A.26) A.6.2 Pseudo code for iCKAN with implicit time integration and helper network Algorithm 1 iCKAN with implicit time integration and helper network 0: Δt t_t, tF_t, f 1: Retrieve t−1C_t-1, i,t−1U_i,t-1 from history 2: Get trial value of evolution equation i,ttrialD_i,t^trial 3: t←tTtC_t _t^TF_t 4: ¯e,ttrial←i,t−1−1ti,t−1 C_e,t^trial _i,t-1^-1C_tU_i,t-1 5: ψ←cKAN(I1¯e,ttrial,I2¯e,ttrial,I3¯e,ttrial)ψ (I_1 C_e,t^trial,I_2 C_e,t^trial,I_3 C_e,t^trial) KAN for elastic potential 6: ¯trial←2¯e,ttrial∂ψ∂¯e,ttrial ^trial← 2 C_e,t^trial ∂ψ∂ C_e,t^trial 7: ω←cKAN(I1¯trial,J2¯trial,J3¯trial3)ω (I_1 ^trial, J_2 ^trial, [3]J_3 ^trial) KAN for inelastic potential 8: itrial←∂ω∂¯trialD_i^trial← ∂ω∂ ^trial 9: Predict i,t←fNN(i,t−1,t,i,ttrial)L_i,t← f_N(U_i,t-1,C_t,D_i,t^trial) 10: i,t←iiTU_i,t _iL_i^T 11: Get predicted value of evolution equation i,tD_i,t 12: i←∂ω∂¯D_i← ∂ω∂ 13: Get residual of evolution equation: :=i,t−1−i,texp(−2Δti,t)i,t=!r:=C_i,t-1-U_i,t (-2\, t\,D_i,t)U_i,t !=0 14: ¯e,t←i,t−1ti,t C_e,t _i,t^-1C_tU_i,t 15: Compute invariants: I1,e,I2,e,I3,e←Invariants(¯e,t)I_1,e,I_2,e,I_3,e ( C_e,t) 16: ψ←cKAN(I1,e,I2,e,I3,e)ψ (I_1,e,I_2,e,I_3,e) KAN for elastic potential 17: t←2i−1∂ψ∂¯e,ti−1S_t← 2U_i^-1 ∂ψ∂ C_e,tU_i^-1 18: t←ttP_t _tS_t 19: Update history: store tC_t, i,tU_i,t 20: return tP_t A.7 Candidate functions for symbolification Candidate functions that are convex and non-decreasing on the interval [0,∞)[0,∞) [Thakolkaran et al., 2025]: x, x1.5x^1.5, x2x^2, x2.5x^2.5, x3x^3, x3.5x^3.5, x4x^4, x4.5x^4.5, x5x^5 exe^x, log(ex) (e^x), log(ex)2 (e^x)^2, log(ex)3 (e^x)^3 Candidate functions that are convex, zero-valued and non-negative on the interval (−∞,∞)(-∞,∞): |x||x|, |x|1.5|x|^1.5, x2x^2, |x|2.5|x|^2.5, |x|3|x|^3, |x|3.5|x|^3.5, x4x^4, |x|4.5|x|^4.5, |x|5|x|^5 ex2−1e^x^2-1, ex4−1e^x^4-1 cosh(x)−1cosh(x)-1, cosh(x2)−1cosh(x^2)-1 log(cosh(x)) (cosh(x)), log(cosh(x2)) (cosh(x^2)) A.8 Additional results of iCKANs on synthetic data In this section, we show the additional results of the iCKANs trained on synthetic data. Using the explicit time integration, we examine the options where the input arguments of the KAN for inelastic potential is augmented with the negative stress invariants or additional pricipal invariants. We further examine the performance of iCKANs with implicit time intergration scheme Appendix A.6. We provide the details on the performance of each variant and compare the convergence between the explicit and implicit iCKANs. Explicit iCKAN: Full result. The prediction of the symbolified version of the explicit iCKAN on the entire dataset is shown in Figure A.2 tload=0.2t_load=0.2 stload=0.5t_load=0.5 stload=1t_load=1 s Figure A.2: Predicted first Piola-Kirchhoff stress by the symbolified explicit iCKAN. The gray area shows the range of the training data and the dots with corresponding color indicate the reference data. Explicit iCKAN: Adding negative stress invariants to arguments of inelasitc potential. We could enhance the arguments for the inelastic potential with the negative stress invariants, as shown in Equation A.11. This approach provides higher expressibility of the network, but unfortunately shows lower stability in prediction. To achieve more interpretable symbolifc expressions of the inelasitc potential we could symbolify the combined activations for the positive and negative parts of each input argument together, i.e. y1symb=y1+(I^1¯)+y1−(−I^1¯),y2symb=y2+(J^2¯)+y2−(−J^2¯),y3symb=y3+(J^3¯)+y3−(−J^3¯). split&y^symb_1=y_1^+( I_1 )+y_1^-(- I_1 )\,,\\ &y^symb_2=y_2^+( J_2 )+y_2^-(- J_2 )\,,\\ &y^symb_3=y_3^+( J_3 )+y_3^-(- J_3 )\,. split (A.27) The candidate functions are general convex functions that are not necessarily non-decreasing. The architecture of the KAN for inelastic potential before and after the symbolification is shown in Figure A.3. The last three activations are pruned to be zero since their contribution to the output is included in the first three activations. It is to note, that the ℋH-transformation has to be applied to the network outcome to achieve convex, non-negative and zero-valued at origin inelastic potential. overpic[width=136.5746pt]figures/KAN_61.png (-1.0,17.0) [width=20.48514pt]figures/spline/omega6/sp_0_0_0.png (16.0,17.0) [width=20.48514pt]figures/spline/omega6/sp_0_1_0.png (33.0,17.0) [width=20.48514pt]figures/spline/omega6/sp_0_2_0.png (51.0,17.0) [width=20.48514pt]figures/spline/omega6/sp_0_3_0.png (68.0,17.0) [width=20.48514pt]figures/spline/omega6/sp_0_4_0.png (86.0,17.0) [width=20.48514pt]figures/spline/omega6/sp_0_5_0.png (3.0,-10.0)$ I_1 $ (18.0,-10.0)$ J_2 $ (35.0,-10.0)$ J_3 $ (51.0,-10.0)$- I_1 $ (70.0,-10.0)$- J_2 $ (89.0,-10.0)$- J_3 $ (50.0,64.0)$ω^KAN$ overpic (a) ⟶Symb. Symb. overpic[width=136.5746pt]figures/KAN_61_pruned.png (-1.0,17.0) [width=20.48514pt]figures/spline_symb/omega6/sp_0_0_0.png (17.0,17.0) [width=20.48514pt]figures/spline_symb/omega6/sp_0_1_0.png (35.0,17.0) [width=20.48514pt]figures/spline_symb/omega6/sp_0_2_0.png (3.0,-10.0)$ I_1 $ (18.0,-10.0)$ J_2 $ (35.0,-10.0)$ J_3 $ (51.0,-10.0)$- I_1 $ (70.0,-10.0)$- J_2 $ (89.0,-10.0)$- J_3 $ (50.0,64.0)$ω^KAN$ overpic (b) Figure A.3: KAN for inelastic potential using both positive and negative stress invariants as arguments (a) before and (b) after the symbolification of the activations. Explicit iCKAN: Adding princiap invariants to arguments of inelasitc potential. We could enhance the arguments of the KAN for inelastic potential with the principal invariants, as shown in Equation A.14. Implicit iCKAN. We consider the implicit iCKAN variant using three modified stress invariants as inputs to the KAN for inelastic potential (Equation 26). The accuracy of the inelastic stretch predicted by the helper LTC network is compared with the Newton–Raphson solution of the implicit evolution equation in Figure A.4 for the case of F11max=1.3F_11 =1.3 and F˙11=0.6 F_11=0.6 s−1-1. Although trained only on the gray-marked range, the LTC shows good generalization to larger deformation levels. t[s] Figure A.4: Stress prediction of the implicit iCKAN. During training, the evolution equation is solved with a trainable LTC helper network. The gray region indicate the range of training data. Predictions using the trained LTC are compared with solutions obtained by a Newton–Raphson solver. Comparison between the explicit and implicit iCKAN. For the case where three modified stress invariants are used as input arguments for the KAN for inelastic potential, the convergence behavior of the explicit and implicit time integration schemes is compared in Figure A.5. (a) Explicit training(b) Implicit training with helper LTC Figure A.5: Training loss of the iCKAN model with L1 regularization using (a) the explicit integration scheme and (b) the implicit integration scheme with helper LTC to solve the evolution equation. Detailed results. Table A.1 shows the detailed results for both explicit and implicit iCKAN before and after symbolification. The quantities ℒstressL_stress and ℒtestL_test denote the NMSE of the predicted stress on the training and full datasets, respectively. ℒsymbL_symb is the NMSE of the symbolified network on the full dataset, and ℒevoL_evo is the loss of the implicit evolution equation when a helper LTC is used during training. Table A.1: Detailed results of training and testing of iCKANs on the synthetic dataset. iCKAN ωKAN(∙)ω^KAN( ) ℒstressL_stress ℒevoL_evo ℒtestL_test ℒsymbL_symb Explicit (I^1¯,J^2¯,J^3¯)( I_1 , J_2 , J_3 ) 2.0⋅10−52.0· 10^-5 - 6.3⋅10−46.3· 10^-4 3.8⋅10−43.8· 10^-4 Explicit (I^1¯,J^2¯,J^3¯,I^2¯,I^3¯)( I_1 , J_2 , J_3 , I_2 , I_3 ) 4.2⋅10−44.2· 10^-4 - 9.7⋅10−49.7· 10^-4 - Implicit (LTC) (I^1¯,J^2¯,J^3¯)( I_1 , J_2 , J_3 ) 6.9⋅10−56.9· 10^-5 2.9⋅10−82.9· 10^-8 7.3⋅10−47.3· 10^-4 - Implicit (NR) (I^1¯,J^2¯,J^3¯)( I_1 , J_2 , J_3 ) - - 3.4⋅10−43.4· 10^-4 3.9⋅10−43.9· 10^-4 A.9 Additional data for training iCKANs on experimental data Hyperparameters. Table A.2 lists the hyperparameters of the iCKANs from Section 5 for the material model discovery of VHB 4910 and VHB 4905. Hyperparameter Value VHB 4910 VHB 4905 elastic potential inelastic potential elastic potential inelastic potential Topology [3,2,1] [3,2,1] [4,2,1] [3,2,1] Order of splines 3 3 3 3 Grid intervals 1 1 1 1 Optimizer AMSGrad AMSGrad Scheduler Cyclic LR Cyclic LR Base Learning rate 5⋅10−35· 10^-3 1⋅10−41· 10^-4 Max Learning rate 5⋅10−25· 10^-2 1⋅10−31· 10^-3 Clip gradient norm 0.10.1 0.10.1 L1 regularization magnitude 1⋅10−51· 10^-5 1⋅10−41· 10^-4 Table A.2: Hyperparameters for the iCKAN models used in the numerical examples in Section 5 for VHB 4910 and VHB 4905. Detailed results. Table A.3 shows the detailed results of iCKANs performance on VHB 4910 and VHB 4905 polymer. The stress loss before symbolification is denoted as ℒL and the stress loss after symbolification is denoted as ℒ L. Figure A.6 shows the loss convergence of training iCKANs from Section 5 for the material model discovery of VHB 4910 and VHB 4905. Table A.3: Detailed results of training and testing of iCKANs on VHB 4910 and VHB 4905 dataset. ℒtrainL_train ℒtestL_test ℒ^train L_train ℒ^test L_test VHB 4910 7.6⋅10−57.6· 10^-5 1.3⋅10−31.3· 10^-3 1.9⋅10−41.9· 10^-4 1.9⋅10−31.9· 10^-3 VHB 4905 2.2⋅10−42.2· 10^-4 6.5⋅10−46.5· 10^-4 4.1⋅10−44.1· 10^-4 6.9⋅10−46.9· 10^-4 Loss (a) VHB 4910 Loss (b) VHB 4905 Figure A.6: Losses during training of the iCKAN for the experimental data of (a) VHB 4910 and (b) VHB4905 polymer. The losses are plotted on a logarithmic scale. For training, 1300 epochs were used. Appendix B Declarations B.1 Acknowledgements This work was supported by the Emmy Noether Grant 533187597 by the Deutsche Forschungsgemeinschaft to Kevin Linka. Christian Cyron greatfully acknowledges support of the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation programme (grant agreement No 101167207 / MechVivo). B.2 Conflict of interest The authors of this work certify that they have no affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; participation in speakers’ bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. B.3 Code availability The source code will be made publicly available upon acceptance of the manuscript. B.4 Contributions by the authors Chenyi Ji: Writing – original draft, Writing – review & editing, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Kian P. Abdolazizi: Writing – review & editing, Supervision, Software, Methodology, Data curation, Conceptualization. Hagen Holthusen: Writing – review & editing, Supervision, Methodology, Conceptualization. Chrstian J. Cyron: Writing – review & editing, Methodology, Conceptualization. Kevin Linka: Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization. B.5 Declaration of generative AI and AI-assisted technologies in the manuscript preparation process During the preparation of this work the authors used OpenAI’s ChatGPT in order to refine the language. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article. References K. P. Abdolazizi, K. Linka, and C. J. Cyron (2024) Viscoelastic constitutive artificial neural networks (vcanns)–a framework for data-driven anisotropic nonlinear finite viscoelasticity. Journal of computational physics 499, p. 112704. Cited by: §5.2. K. P. Abdolazizi, R. C. Aydin, C. J. Cyron, and K. Linka (2025) Constitutive Kolmogorov–Arnold Networks (CKANs): Combining accuracy and interpretability in data-driven material modeling. Journal of the Mechanics and Physics of Solids 203, p. 106212. External Links: Document, 2502.05682, ISSN 00225096, Link Cited by: §A.2, §A.4, §1, §3.1, §3.2.1, §3.2, §3, §4, §5.2.2, §6, §6, §6, §6, §6. R. Abdusalamov, M. Hillgärtner, and M. Itskov (2023) Automatic generation of interpretable hyperelastic material models by symbolic regression. International Journal for Numerical Methods in Engineering 124 (9), p. 2093–2104. Cited by: §1. D. W. Abueidda, P. Pantidis, and M. E. Mobasher (2025) Deepokan: deep operator network based on kolmogorov arnold networks for mechanics problems. Computer Methods in Applied Mechanics and Engineering 436, p. 117699. Cited by: §1. B. Amos, L. Xu, and J. Z. Kolter (2017) Input convex neural networks. In International conference on machine learning, p. 146–155. Cited by: §3.1.2. F. As’ ad and C. Farhat (2023) A mechanics-informed neural network framework for data-driven nonlinear viscoelasticity. In AIAA SCITECH 2023 Forum, p. 0949. Cited by: §A.6.1, §1. C. Benoit (1924) Note sur une méthode de résolution des équations normales provenant de l’application de la méthode des moindres carrés à un système d’équations linéaires en nombre inférieur à celui des inconnues (procédé du commandant cholesky). Bulletin géodésique 2 (1), p. 67–77. Cited by: §A.6.1. B. Boes, J. Simon, and H. Holthusen (2026) Accounting for plasticity: an extension of inelastic constitutive artificial neural networks. European Journal of Mechanics - A/Solids 117, p. 105998. External Links: ISSN 0997-7538 Cited by: §1, §6. G. Bomarito, T. Townsend, K. Stewart, K. Esham, J. Emery, and J. Hochhalter (2021) Development of interpretable, data-driven plasticity models with symbolic regression. Computers & Structures 252, p. 106557. Cited by: §1. F. Dammaß, K. A. Kalina, and M. Kästner (2025) Neural networks meet phase-field: a hybrid fracture model. Computer Methods in Applied Mechanics and Engineering 440, p. 117937. Cited by: §1. T. Deschatre and X. Warin (2025) Input convex kolmogorov arnold networks. arXiv preprint arXiv:2505.21208. Cited by: §3.1.2. M. Fernández, F. Fritzen, and O. Weeger (2022) Material modeling for parametric, anisotropic finite strain hyperelasticity based on machine learning with application in optimization of metamaterials. International Journal for Numerical Methods in Engineering 123 (2), p. 577–609. Cited by: §1. M. Fernández, M. Jamshidian, T. Böhlke, K. Kersting, and O. Weeger (2021) Anisotropic hyperelastic constitutive models for finite deformations combining material theory and data-driven approaches with application to cubic lattice metamaterials. Computational Mechanics 67 (2), p. 653–677. Cited by: §1. M. Flaschel, S. Kumar, and L. De Lorenzis (2021) Unsupervised discovery of interpretable hyperelastic constitutive laws. Computer Methods in Applied Mechanics and Engineering 381, p. 113852. Cited by: §1. M. Flaschel, S. Kumar, and L. De Lorenzis (2023) Automated discovery of generalized standard material models with euclid. Computer Methods in Applied Mechanics and Engineering 405, p. 115867. Cited by: §1, §6. M. Flaschel, P. Steinmann, L. De Lorenzis, and E. Kuhl (2025) Convex neural networks learn generalized standard material models. Journal of the Mechanics and Physics of Solids 200, p. 106103. Cited by: §1. M. Flaschel (2023) Automated discovery of material models in continuum solid mechanics. Ph.D. Thesis, ETH Zurich. Cited by: §1. P. Flory (1961) Thermodynamic relations for high elastic materials. Transactions of the Faraday Society 57, p. 829–838. Cited by: §3.2.1. J. N. Fuhg, G. A. Padmanabha, N. Bouklas, B. Bahmani, W. Sun, N. N. Vlassis, M. Flaschel, P. Carrara, and L. De Lorenzis (2024) A review on data-driven constitutive laws for solids. arXiv preprint arXiv:2405.03658. Cited by: §1. P. Germain, Q. S. Nguyen, and P. Suquet (1983) Continuum thermodynamics. Journal of applied mechanics 50, p. 1010–1020. Cited by: §2. B. Guo, Z. Lin, and Q. He (2025) History-aware neural operator: robust data-driven constitutive modeling of path-dependent materials. arXiv preprint arXiv:2506.10352. Cited by: §1. B. Halphen and Q. S. Nguyen (1975) Sur les matériaux standard généralisés. Journal de mécanique 14 (1), p. 39–63. Cited by: §1, §2. S. Hartmann and P. Neff (2003) Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near-incompressibility. International journal of solids and structures 40 (11), p. 2767–2791. Cited by: §3.2.1. R. Hasani, M. Lechner, A. Amini, D. Rus, and R. Grosu (2021) Liquid time-constant networks. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 35, p. 7657–7666. Cited by: §A.6.1. H. Holthusen, T. Brepols, K. Linka, and E. Kuhl (2025) Automated model discovery for tensional homeostasis: constitutive machine learning in growth and remodeling. Computers in biology and medicine 186, p. 109691. Cited by: §1, §3.2.4, §6. H. Holthusen and E. Kuhl (2026) A complement to neural networks for anisotropic inelasticity at finite strains. Computer Methods in Applied Mechanics and Engineering 450, p. 118612. External Links: ISSN 0045-7825 Cited by: §A.6.1, §A.6, §1, §3.2.4, §6. H. Holthusen, L. Lamm, T. Brepols, S. Reese, and E. Kuhl (2024a) Polyconvex inelastic constitutive artificial neural networks. PAMM 24 (3), p. e202400032. Cited by: §3.2.1. H. Holthusen, L. Lamm, T. Brepols, S. Reese, and E. Kuhl (2024b) Theory and implementation of inelastic constitutive artificial neural networks. Computer Methods in Applied Mechanics and Engineering 428, p. 117063. Cited by: §2, §2, §2, §3.2.1. H. Holthusen, K. Linka, E. Kuhl, and T. Brepols (2026) A generalized dual potential for inelastic constitutive artificial neural networks: a jax implementation at finite strains. Journal of the Mechanics and Physics of Solids 206, p. 106337. External Links: ISSN 0022-5096 Cited by: §A.3, §A.5, §A.5, §1, §2, §2, §2, §3.2.2, §3.2.4, §3.2.4, §5.2, §6, §6, §6, §6. H. Holthusen, C. Rothkranz, L. Lamm, T. Brepols, and S. Reese (2023) Inelastic material formulations based on a co-rotated intermediate configuration—application to bioengineered tissues. Journal of the Mechanics and Physics of Solids 172, p. 105174. External Links: ISSN 0022-5096 Cited by: Figure 2, §2, §2. G. A. Holzapfel (2000) Nonlinear solid mechanics: a continuum approach for engineering science. John Wiley & Sons, Chichester. External Links: ISBN 0-471-82319-8 Cited by: §A.2. K. Hornik, M. Stinchcombe, and H. White (1989) Multilayer feedforward networks are universal approximators. Neural networks 2 (5), p. 359–366. Cited by: §1. T. Hospedales, A. Antoniou, P. Micaelli, and A. Storkey (2021) Meta-learning in neural networks: a survey. IEEE transactions on pattern analysis and machine intelligence 44 (9), p. 5149–5169. Cited by: §6. M. Hossain, D. K. Vu, and P. Steinmann (2012) Experimental study and numerical modelling of vhb 4910 polymer. Computational Materials Science 59, p. 65–74. Cited by: §5.2.1, §5. B. Hudobivnik and J. Korelc (2016) Closed-form representation of matrix functions in the formulation of nonlinear material models. Finite Elements in Analysis and Design 111, p. 19–32. Cited by: §A.5, §3.2.4. A. A. Jadoon, K. A. Meyer, and J. N. Fuhg (2025) Automated model discovery of finite strain elastoplasticity from uniaxial experiments. Computer Methods in Applied Mechanics and Engineering 435, p. 117653. Cited by: §1. T. Ji, Y. Hou, and D. Zhang (2024) A comprehensive survey on kolmogorov arnold networks (kan). arXiv preprint arXiv:2407.11075. Cited by: §1. R. E. Jones, A. L. Frankel, and K. Johnson (2022) A neural ordinary differential equation framework for modeling inelastic stress response via internal state variables. Journal of Machine Learning for Modeling and Computing 3 (3). Cited by: §1. E. Kabliman, A. H. Kolody, J. Kronsteiner, M. Kommenda, and G. Kronberger (2021) Application of symbolic regression for constitutive modeling of plastic deformation. Applications in Engineering Science 6, p. 100052. Cited by: §1. K. A. Kalina, J. Brummund, W. Sun, and M. Kästner (2025) Neural networks meet anisotropic hyperelasticity: a framework based on generalized structure tensors and isotropic tensor functions. Computer Methods in Applied Mechanics and Engineering 437, p. 117725. Cited by: §6. T. Kirchdoerfer and M. Ortiz (2016) Data-driven computational mechanics. Computer Methods in Applied Mechanics and Engineering 304, p. 81–101. Cited by: §1. E. Kiyani, K. Shukla, J. F. Urbán, J. Darbon, and G. E. Karniadakis (2025) Optimizing the optimizer for physics-informed neural networks and kolmogorov-arnold networks. Computer Methods in Applied Mechanics and Engineering 446, p. 118308. Cited by: §1. A. N. Kolmogorov (1961) On the representation of continuous functions of several variables by superpositions of continuous functions of a smaller number of variables. American Mathematical Society. Cited by: §3.1. J. Li, Z. Guan, J. Chen, and H. Jin (2025) A long short-term memory-based constitutive modeling framework for capturing strain path dependence in plastic deformation. Mechanics of Materials 205, p. 105325. Cited by: §1. Z. Liao, M. Hossain, X. Yao, M. Mehnert, and P. Steinmann (2020) On thermo-viscoelastic experimental characterization and numerical modelling of vhb polymer. International Journal of Non-Linear Mechanics 118, p. 103263. Cited by: §5.2.2, §5. L. Linden, D. K. Klein, K. A. Kalina, J. Brummund, O. Weeger, and M. Kästner (2023) Neural networks meet hyperelasticity: a guide to enforcing physics. Journal of the Mechanics and Physics of Solids 179, p. 105363. Cited by: §1. K. Linka, M. Hillgärtner, K. P. Abdolazizi, R. C. Aydin, M. Itskov, and C. J. Cyron (2021) Constitutive artificial neural networks: a fast and general approach to predictive data-driven constitutive modeling by deep learning. Journal of Computational Physics 429, p. 110010. Cited by: §1. K. Linka and E. Kuhl (2023) A new family of constitutive artificial neural networks towards automated model discovery. Computer Methods in Applied Mechanics and Engineering 403, p. 115731. Cited by: §1. Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljačić, T. Y. Hou, and M. Tegmark (2024) Kan: kolmogorov-arnold networks. arXiv preprint arXiv:2404.19756. Cited by: §1, §3, §4. F. Masi, I. Stefanou, P. Vannucci, and V. Maffi-Berthier (2021) Thermodynamics-based artificial neural networks for constitutive modeling. Journal of the Mechanics and Physics of Solids 147, p. 104277. Cited by: §1. M. Maurizi, C. Gao, and F. Berto (2022) Predicting stress, strain and deformation fields in materials and structures with graph neural networks. Scientific reports 12 (1), p. 21834. Cited by: §1. V. K. Narouie, J. Urrea-Quintero, F. Cirak, and H. Wessels (2026) Unsupervised constitutive model discovery from sparse and noisy data. Computer Methods in Applied Mechanics and Engineering 452, p. 118722. Cited by: §1. A. Polo-Molina, D. Alfaya, and J. Portela (2024) MonoKAN: certified monotonic kolmogorov-arnold network. arXiv preprint arXiv:2409.11078. Cited by: §3.1.1, §3.2.1. M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, p. 686–707. Cited by: §1. S. J. Reddi, S. Kale, and S. Kumar (2019) On the convergence of adam and beyond. arXiv preprint arXiv:1904.09237. Cited by: §5. M. Rosenkranz, K. A. Kalina, J. Brummund, W. Sun, and M. Kästner (2024) Viscoelasticty with physics-augmented neural networks: model formulation and training methods without prescribed internal variables. Computational Mechanics 74 (6), p. 1279–1301. Cited by: §A.6.1, §1, §1. L. N. Smith (2017) Cyclical learning rates for training neural networks. In 2017 IEEE winter conference on applications of computer vision (WACV), p. 464–472. Cited by: §5. V. Tac, F. S. Costabal, and A. B. Tepole (2022) Data-driven tissue mechanics with polyconvex neural ordinary differential equations. Computer Methods in Applied Mechanics and Engineering 398, p. 115248. Cited by: §1. V. Taç, M. K. Rausch, F. S. Costabal, and A. B. Tepole (2023) Data-driven anisotropic finite viscoelasticity using neural ordinary differential equations. Computer methods in applied mechanics and engineering 411, p. 116046. Cited by: §1. M. Tacke, M. Busch, K. Bali, K. Abdolazizi, K. Linka, C. Cyron, and R. Aydin (2025) Constitutive scientific generative agent (csga): leveraging large language models for automated constitutive model discovery. Machine Learning for Computational Science and Engineering 1 (1), p. 23. Cited by: §6. P. Thakolkaran, Y. Guo, S. Saini, M. Peirlinck, B. Alheit, and S. Kumar (2025) Can kan cans? input-convex kolmogorov-arnold networks (kans) as hyperelastic constitutive artificial neural networks (cans). Computer Methods in Applied Mechanics and Engineering 443, p. 118089. Cited by: §A.2, §A.4, §A.4, §A.7, §1, §3.1.1, §3.1.1, §3, §6, §6. D. Versino, A. Tonda, and C. A. Bronkhorst (2017) Data driven modeling of plastic deformation. Computer Methods in Applied Mechanics and Engineering 318, p. 981–1004. Cited by: §1. Y. Wang, J. Sun, J. Bai, C. Anitescu, M. S. Eshaghi, X. Zhuang, T. Rabczuk, and Y. Liu (2025) Kolmogorov–arnold-informed neural network: a physics-informed deep learning framework for solving forward and inverse problems based on kolmogorov–arnold networks. Computer Methods in Applied Mechanics and Engineering 433, p. 117518. Cited by: §1. H. You, Q. Zhang, C. J. Ross, C. Lee, M. Hsu, and Y. Yu (2022) A physics-guided neural operator learning approach to model biological tissues from digital image correlation measurements. Journal of Biomechanical Engineering 144 (12), p. 121012. Cited by: §1.