Elucidating the auxetic behavior of cementitious cellular Elucidating the auxetic behavior of cementitious cellular composites using finite element analysis and interpretable composites using finite element analysis and interpretable machine learning machine learning

FEA)-based framework with ML to predict the Poisson’s ratios. In particular, the FEA-based approach is used to generate a dataset containing 850 combinations of different mesoscale architectural void features. The dataset is leveraged to develop an ML-based prediction tool using a feed-forward multilayer perceptron-based neural network (NN) approach which shows excellent prediction efﬁcacy. To shed light on the relative inﬂuence of the design parameters on the auxetic behavior of the ACCCs, Shapley additive explanations (SHAP) is employed, which establishes the volume fraction of voids as the most inﬂuential parameter in inducing auxetic behavior. Overall, this paper develops an efﬁcient approach to evaluate geometry-dependent auxetic behaviors for cementitious materials which can be used as a starting point toward the design and development of auxetic behavior in cemen-titious composites


a b s t r a c t
With the advent of 3D printing, auxetic cellular cementitious composites (ACCCs) have recently garnered significant attention owing to their unique mechanical performance. To enable seamless performance prediction of the ACCCs, interpretable machine learning (ML)-based approaches can provide efficient means. However, the prediction of Poisson's ratio using such ML approaches requires large and consistent datasets which is not readily available for ACCCs. To address this challenge, this paper synergistically integrates a finite element analysis (FEA)-based framework with ML to predict the Poisson's ratios. In particular, the FEA-based approach is used to generate a dataset containing 850 combinations of different mesoscale architectural void features. The dataset is leveraged to develop an ML-based prediction tool using a feed-forward multilayer perceptron-based neural network (NN) approach which shows excellent prediction efficacy. To shed light on the relative influence of the design parameters on the auxetic behavior of the ACCCs, Shapley additive explanations (SHAP) is employed, which establishes the volume fraction of voids as the most influential parameter in inducing auxetic behavior. Overall, this paper develops an efficient approach to evaluate geometry-dependent auxetic behaviors for cementitious materials which can be used as a starting point toward the design and development of auxetic behavior in cementitious composites.

Introduction
Ordinary Portland cement (OPC)-based cementitious composites are one of the most widely used building materials in the world [1]. However, the performance of concrete falls short under tension, fracture, and high strain rate conditions [2][3][4][5][6]. Improvement of concrete properties through tuning of micro and mesostructure has garnered significant attention over the previous several decades. These research efforts have focused on performance-enhancement of concrete through optimizing the pore structure [7][8][9][10], modifications in the aggregate type [11][12][13], tuning of packing density [14][15][16], and incorporating new reinforcing phases in concrete [17][18][19][20][21]. With the advent of 3D printing technology, it is now possible to achieve performance enhancement by tuning the geometrical architecture in the mesoscale. In particular, auxetic cementitious composites with negative Poisson's ratios have been in the spotlight recently owing to their excellent mechanical performance [22][23][24]. Implementation of axial tension in a material with auxetic behavior leads to an increase in the cross-section, while axial compression induces a transversal contraction [25]. Auxetic structures have been shown to offer significant energy absorption under dynamic loading conditions due to such unique mechanical behavior [26][27][28][29][30][31]. Auxetic behavior has been studied in detail for 3D printed polymers [32][33][34][35], metals [36][37][38], and ceramics [39,40]. There are different approaches to achieve auxetic behavior such as implementation of re-entrant honeycombs [25,33,41], rigid body rotation structures [42][43][44], crumpled or perforated sheet structures [45][46][47], chiral structures [48][49][50], folded structures [51][52][53], double arrowed structures [26,27,[29][30][31] or buckling-induced structures [37,54,55]. The mechanical response of a three-dimensional double arrowed (or double -V) auxetic structure under axial loading has been studied in [56] where the theoretical and numerical quasistatic analyses are carried out and validated with the experimental results with good accuracy. In addition, the geometrically nonlinear mechanical properties of the auxetic double-V structure are also extended to large strains using a large beam deflection model where good agreement was found between the simulated and experimental results [57]. In the study by Lin et al. [58], chiral mechanical metamaterials have been found to exhibit a unique compression-twisting coupling effect under external compressive loading. The effects of chirality on the buckling strength and buckling mode were also evaluated in the study. These chiral metamaterials have been shown to intrinsically undergo self-twist and control their global rotation [59]. The effect of global rotation is observed both for small systems such as single structural units and for complex structures. These chiral metamaterials owing to their compression/thermal-torsion effects demonstrated not only negative Poisson's ratio (PR) but also the negative coefficient of thermal expansion [60]. Thus, the coupling design method for metamaterials provides a good design strategy for a wide range of engineering applications. Besides, because of the unusual but desirable properties such as enhanced indentation resistance, variable permeability, superior shear resistance, high fracture toughness, and damping and sound absorption [61,62], auxetic materials are used for a variety of applications. These include but are not limited to sports and safety equipment, shock and soundabsorbing materials for vehicles and aircraft, seals, filters, prostheses, and textiles [25,61,62].
While 3D printed cellular polymers and metals showed impressive auxetic performance with high fracture toughness [33,63,64] and impact resistance [65][66][67][68][69][70][71], the brittleness and low deformability of cellular cementitious composites in comparison with materials like polymers make it challenging to achieve auxetic behavior. In fact, the implementation of auxetic structures by using strategically placed holes of different shapes in concretes and cementitious materials is relatively new. Thanks to the significant improvements in the 3D printing technology for cementitious materials that have made this conceptual idea a reality. Auxetic cementitious composites have the potential to offer higher energy dissipation and higher impact resistance to build resilient civil infrastructure. For example, having higher energy dissipation in concretes can help toward the design of efficient earthquakeresistant structures. So far, auxetic concretes and cementitious materials have not been implemented in civil infrastructure, and only a limited number of experimental as well as simulation studies exist. To induce auxetic behavior in cementitious composites, Xu et al. recently implemented periodic unit cells of cellular cementitious composites with elliptical voids [22,23] and the mechanical behavior was evaluated using both experiments and simulations. Four-point-bending, uniaxial compression, and cyclic loading tests on cube and beam specimens, respectively were performed to confirm the simulation results. The digital image correlation (DIC) method was used to determine the transverse and the axial strain in order to identify the Poisson's ratio and with it the auxetic behavior [22]. With the same structure type of alternate ellipses, different geometries with different unit cell sizes and ratios between the major and the minor axis of the ellipse were also investigated.
While these previous works [22][23][24] show that auxetic behavior of cellular cementitious composites is generally possible and explain all the benefits and potentials, implementation of this idea on a large industrial scale has not been realized yet in the civil infrastructure primarily due to the lack of reliable performance standards and unavailability of a robust yet user-friendly predictive software. Thus, the development of a simple and yet robust as well as user-friendly design tool would provide enormous value to the designers and decision-makers as they would be able to easily leverage such tools for performance tuning and risk assessment which is challenging to do with computationally exhaustive and complex finite element analysis (FEA)-based simulation tools developed so far. Hence this paper is aimed at efficient interpretable machine learning (ML)-based performance prediction for such cellular cementitious composites. Such ML-based predictive algorithms can be packaged into user-friendly software that can be easily leveraged by engineers and decision-makers. Besides, the implementation of such an interpretable ML-based approach can open up possibilities to uncover various unexplored domains of design strategies for these auxetic cellular cementitious composites. However, the development of a robust ML-based prediction tool requires a large, consistent, and representative dataset that is currently not available for auxetic cellular cementitious composites. Moreover, owing to the highly complex nature (which is also referred to as black-box) of ML techniques such as neural network (NN) [72], its indiscriminate use can yield non-physical solutions [73,74]. Therefore, to address the aforementioned challenges, this study implements high-throughput FE simulations to obtain a large dataset of Poisson's ratio for varying mesoscale architectural features. The large dataset is leveraged for ML-based performance prediction. Integration of ML with FEA-based numerical simulation helps to maintain the fundamental laws of physics and to avoid nonphysical solutions. This paper implements a feed-forward multilayer perceptron-based NN approach primarily because of its potential high accuracy of prediction as demonstrated by several materials-design studies [75][76][77][78][79][80][81][82]. Additionally, this paper also implements SHapley Additive exPlanations (SHAP) [83] to interpret the results obtained from NN and to elucidate the relative sensitivity of the design parameters towards the predicted Poisson's ratio. Overall, the synergistically integrated high throughput FEA-based simulations and ML-based approach, presented in this paper, is expected to enable materials engineers and decisionmakers to make intelligent, informed decisions regarding the selection of mesoscale architectural features in cellular cementitious composites to meet the desired performance goals. Thus, the design tools presented in this paper can accelerate the acceptance and utilization of auxetic cementitious cellular composites for a wide range of structural applications.

Response prediction using finite element analysis
This section explains how the effective Poisson's ratio of the cellular cementitious composites is computed using finite element analysis. The auxetic behavior in these cellular cementitious composites is characterized by the negative Poisson's ratio. The forthcoming sub-sections first detail the numerical simulation methodology followed by response evaluations and large dataset generation.

Simulation methodology
In this section, the numerical simulation framework towards the prediction of Poisson's ratio of the cellular cementitious composites with variating shape and dimension of voids is elaborated. The framework is schematically illustrated in Fig. 1.
The first step involves generating the periodic geometry for representative unit cells with elliptical voids. Following the geometry generation, the unit cells are meshed. Mesh-convergence studies (please refer to the supplementary material) are performed to determine the optimal mesh for converged solutions. Afterward, the Dirichlet boundary condition is chosen to mimic the free edges uniaxial compression experiment. A uniaxial strain at a finite rate is applied, and longitudinal strain and transverse strain responses are obtained. Using a homogenization procedure and the linear relationship between the longitudinal strain and the transverse strain, Poisson's ratio is then computed. The analysis is carried out in ABAQUS TM using python scripts for geometry generation and followed by a post-processing module coded in MATLAB. The forthcoming sub-sections elaborate on each component of the simulation framework.

Geometry generation approach, model geometry, and parametric variations
The unit cell is generated by first selecting the desired number of voids along the X-axis and the Y-axis. The center points of the voids are placed at equal spacings based on the size of the unit cell. Once, the points have been generated, a unit circle is filled in with a given diameter. For an aspect ratio greater than 1, the ellipse shape is generated by varying the unit circle with scale factors that correspond to the major and minor axes length or the major axis length and the aspect ratio. As in this study, the elliptical shapes are arranged alternatively (as shown in Fig. 2), therefore, two angles are considered at 0 degrees and 90 degrees. These angles are necessary to rotate the voids accordingly. Since, the aspect ratio considered in this study is greater than or equal to 1, the maximum diameter (or the maximum major axis length) is varied until the minimum specified spacing between voids is achieved (which is 2 mm, in this study). Fig. 2 shows one representative geometry with the distribution of voids along the x and y-axis. Please refer to the supplementary material for additional geometries with varying input features. Table 1 describes the range of the parametric variations considered in this study. The ranges have been selected judiciously with a maximum void content of 50%. The aspect ratio, which is the ratio of the major axis to the minor axis is taken between 1 and 5. The aspect ratio adopted in the literature ranges from 1 to 3 [22,23]. The minimum length of the major axis is 6 mm whereas the maximum length of the major axis is 38.3 mm. Considering the median  diameter of 0.6 mm for sand in mortar, the minimum void length of 6 mm was chosen as smaller voids are challenging to create using standard 3D printers. The minimum spacing of 2 mm is adopted between the voids. The maximum number of voids along the X-axis and Y-axis considered in this study is 8. The same value has also been adopted in the previous literature [22,23]. The void geometry descriptors were carefully chosen considering the manufacturability aspect using 3D printing technology. Similar cementitious structures were successfully printed by Xu et al. [22,23] in their experimental study.

Material properties and boundary conditions
In this study, a linear elastic isotropic material is considered for the mortar with Young's modulus of 25 GPa and Poisson's ratio of 0.18 [84]. This Poisson's ratio is not to be confused with the Poisson's ratio of the auxetic mortars with void content which is to be determined by the simulations. Since the objective of this study is to evaluate the prospects of the negative Poisson's ratio by judiciously placing the elliptical voids in the geometry, the post-peak response and the effect of plasticity are not considered herein. In order to mimic the experimental conditions [22], a Dirichlet boundary condition is imposed in the auxetic structure. While the bottom layer is fixed in all directions, the top layer is fixed in all directions except along the loading direction. The loading at the top part is applied linearly with a prescribed displacement such that it replicates quasi-static loading condition. The influence of the boundary condition for such structure has been addressed in [23]. A similar boundary condition is also considered in the previous literature [22,23] to mimic the experimental uniaxial compression loading.

Effective response prediction
After the generated geometries are successfully meshed and the boundary conditions are applied, the geometries are subjected to uniaxial compression loading at the top edge with the loading speed lying in the quasi-static range (1e-4/s). At every loading step, the longitudinal strain (e L ) and the transverse strain (e T ) is com-puted. From there the effective Poisson's ratio (#) is determined, which is expressed as It needs to be noted that the effective Poisson's ratio is calculated over a selected region of interest (40 mm Â 50 mm) away from the top and bottom boundaries as shown in Fig. 2 using a black dashed rectangle to avoid boundary effects. A similar approach was used for the calculation of the Poisson's ratio in soft materials with 2D digital image correlation [85]. Furthermore, the same concept is applied while calculating the Poisson's ratio of the material using finite element analysis [86].

Simulation results and dataset generation for ML
This section shows the stress-strain results of the finite element analysis. From the illustrations, it is possible to deduce whether auxetic behavior has occurred. After that, the influence of the aspect ratio and the void contents are discussed as an interpretation of the finite element results. Finite element analysis is used later to generate a large and consistent dataset for ML implementation. Fig. 3 shows the stress distribution obtained at different strain states with increasing compressive loading. In voids where the major axis is horizontally aligned, tensile stresses are predominantly present, as indicated by the red color, whereas in voids where the main axis is vertically aligned, stress concentrations are observed along the loading direction, as indicated by the blue color. Similar observations can be made for the strain concentration in Fig. 4. For this particular pattern, high compressive stress and strain is identified for voids closer to the loading and fixed edges, whereas high tensile stress and strain are observed at the center as shown in Fig. 3(c) and 4(c).

Influence of aspect ratio and void contents
To illustrate the influence of the aspect ratio on the auxetic response, Fig. 5 plots the Poisson's ratios with varying aspect ratios and the length of the major axis of the voids for various combinations of the number of voids along the x-axis and y-axis. It is observed that as the number of voids increases, the Poisson's ratio in general decreases and goes into the negative regime. For the cases with a smaller number of voids ( Fig. 5(a)), the distribution is sensitive to the increasing length of the major axis irrespective of the aspect ratio. It can also be detected that irrespective of the aspect ratio, negative Poisson's ratio can be obtained with the increasing length of the major axis. However, it was also observed that the influence of the aspect ratio on Poisson's ratio increases with the increasing number of voids ( Fig. 5(b) and (c)). Thus, with the increasing length of the major axis and an increasing number of voids, the difference of the Poisson's ratios between low and high aspect ratio cases increases significantly. Fig. 6(a), (b), and (c) exhibit the 3D plots of Poisson's ratio distributions with a varying number of voids along the x and y-axis for aspect ratios of 3, 4, and 5 respectively while keeping the length of the major axis fixed at 7.5 mm. It is observed that the negative Poisson's ratio is narrowly distributed towards a higher number of voids (as shown in Fig. 6). Although Poisson's ratio seems to be independent of the aspect ratio for a small number of voids, Poisson's ratio for a high number of voids decreases significantly faster with increasing aspect ratio. For a higher aspect ratio, negative Poisson's ratio can be obtained with the number of voids greater than 4 (as shown in Fig. 6(c)). Using the FE-based simulation framework presented above, a large database of Poisson's ratio val-  ues is generated by varying void contents in terms of the number of voids along the X-axis and Y-axis, aspect ratio, length of the major axis, and volume fraction. This dataset is implemented to develop a machine learning-based performance prediction as explained in the forthcoming section.

Machine learning implementation
3.1. Machine learning techniques and hyper-parameters tuning 3.1.1. Neural network NN is a mathematical model that maps a given set of features, x, to a set of labels (in the case of classification) or a continuous vari-able (in regression), y. This mapping between the features and the output is accomplished through one or more layers of perceptron, which are passed through activation functions to form a feedforward NN. The functional form of an NN can be written as where f N Á ð Þ : R ! R is a continuous bounded function known as the activation function, A i : R d i ! R d iþ1 refers to the transformation matrix which stores the weights in between the two perceptron layers [87]. Due to the universal function approximation capability [88], NNs have been widely adopted to learn the underlying functions associated with a dataset, thereby enabling the discovery of   hidden patterns [89]. This paper uses the 'ReLu' activation function. Further, the weights associated with the perceptrons are obtained as the solution of a constrained optimization given by, where k is the weight associated with the regularization term and g Á ð Þ is the regularization function of the weights, such as L 1 , L 2 , to name a few. The solution can be obtained using any optimizing algorithm, such as stochastic gradient descent or even evolutionary algorithms. Note that in a feed-forward NN, the backpropagation algorithm is leveraged to update the weights. In this study, the Adam algorithm [90] is used as the optimizer.
Due to the non-convexity of this problem, non-unique solutions are obtained for the optimization problem. Moreover, these solutions are highly dependent on the choice of hyperparameters namely, the number of hidden layers and the hidden layer units in each layer. These hyperparameters exhibit significantly higher variations especially when the number of neurons/layers, used in the model, is large. To address these issues, proper regularization techniques should be implemented while training the NN.

Model tuning and cross-validation
To eliminate overfitting, the dataset is split into a training set and a test set. The test set is kept fully hidden while training the model. To assess the performance of the trained NN model, this unseen data (or test set) is used. For this, a k-fold crossvalidation (CV) approach is implemented. . Besides, a nested twolevel CV approach [91] is implemented for hyperparameter tuning. In such an approach, the dataset is first divided into the test set and the training set. Here, 20% of the data is kept hidden from the model as a test set, and the remaining 80% of data in the dataset is used as a training set. While the outer CV is used to compare the performance of the ML algorithm by averaging the values of scores (i.e., R 2 and MSE) obtained from each fold. In order to obtain the appropriate hyperparameters for the NN models, a 5-fold inner cross-validation is adopted for the training set. Such method is successfully implemented in the previous literatures [75,76,92]. A trade-off is sought between accuracy and computation demand while tuning the hyperparameters to obtain the optimized model. This trade-off can be achieved by comparing the performance prediction efficacy for both the training and validation while gradually increasing the number of neurons (degree of complexity). The accuracy of the model is evaluated by calculating the metrics such as linear coefficient of determination (R 2 ) and mean squared error (MSE).

Model interpretability
Complex model techniques such as XGBoost, deep learning, random forest, etc., are often considered as black-box models owing to their low interpretability. Recent technique such as a Shapley additive explanations (SHAP) [83,93] has been developed to interpret these black-box models. In SHAP, the importance of feature 'j' for the output of model f , denoted as / j f ð Þ, is calculated as a weighted sum of the feature's contribution of the model's output f x i ð Þ over all possible feature combinations [94]. This / j f ð Þ is expressed as: Where x j is feature j, S is a subset of features, and p is the number of features in the model. Such approaches have been successfully implemented to interpret the functions learned by the ML models for predicting material properties [75,92,95]. In this study, a SHAP is implemented to interpret the ML results.

Dataset adequacy
To ensure that the dataset used is adequate, there are four factors that the authors have carefully considered in this work. First, the data should be balanced. In this study, the variation of the number of voids, aspect ratio, length of the major axis is incremented consistently such that the variation is uniform and equally important. This is necessary to avoid any imbalance or bias in the dataset. Second, the data should be representative (i.e., that dataset should contain enough information to train the models). During machine learning (ML) training, the data are split into a ratio of 80:20 (training set and test sets). The hyperparameters are tuned using an inner 5-fold cross-validation as defined in section 4.1.2. During each cross-validation fold, the training data is further split into 20% for validation, and 80% for the training set. The average model errors obtained from training error and validation error are represented by taking average values from each fold. The performance of the trained model is evaluated by using the unseen test dataset. If the error between the actual test value and the predicted value is within the acceptable range, the data is considered representative which is true for our current study. Third, the data should be complete. Since the focus of the paper is to predict the Poisson's ratio which is valid till the linear elastic regime, the post-peak response is not implemented. This study covers all the possible ranges for void volume fraction, the number of voids, the aspect ratio, and the length of the major axis for the voids. Fourth, the data should be consistent (i.e., obtained with the same testing protocol). The mechanical analysis for all the parametric variations is modeled using the same procedure where the voids are dispersed as per the number of voids along X-and Y-axes with proper dimensions in the bounding box. The mechanical loading is performed by applying a displacement-controlled uniaxial compressive loading for all the cases. All these simulations are carried in ABAQUS. The same procedure is repeated for all the possible compositions. Thus, the generated data is considered consistent. After considering all the above-mentioned criteria, the number of data points (which is equal to 850 combinations) adopted in the current study is deemed adequate, which can be observed from the performance of the model on the unseen test dataset. To check the performance of each model, the test dataset is not being seen by the model during training and hyper-parameters tuning. The predicted value from the model for the unseen test set exhibits low variation from the actual value as shown later in this paper. This ensures that the data used for training the model are adequate and consistent enough. In general, if the number of data points is increased, the prediction errors are expected to be reduced. However, once the error is saturated, increasing the number of data points does not influence the model performance significantly. Here, in this study, we obtained a good correlation between the predicted and measured responses. Hence, any further increase in the number of data points is not expected to change the model performance significantly. Fig. 7 represents the MSE and R 2 plots with respect to a varying number of neurons for two hidden layers. The number of hidden layers and neurons in NN was based on the hyperparameter optimization in which a trade-off is sought between the minimum MSE and maximum R 2 . From Fig. 7(a), it is observed that that as the number of neurons increases, the MSE value for both training and validation continues to decrease. But, at a neuron equals to 14, the MSE value for the validation set starts to converge, indicating the model has stopped learning and no improvement will be achieved beyond fourteen neurons. Similar behavior is also observed for the R 2 value at neuron equals 14 ( Fig. 7(b)). No significant improvement in the R 2 values is observed with any increase in the number of neurons above 14. Hence, 14 neurons are selected for a trained NN model.

Prediction from trained neural network
Here, the NN model is trained using TensorFlow [96] package in python, where the MSE values between the values obtained from FE simulations and the values predicted by the NN model are minimized to obtain the optimal bias arrays and weight matrices. While training, the parameters (weights and bias) are optimized by implementing back-propagation errors using the Adam algorithm. A learning rate equals 0.001 and 500 epochs are adopted based on the convergence study. The optimized NN features are provided in the supplementary material.
The efficacy of the NN model is evaluated by comparing the predicted values from the trained NN model with the values obtained from FE simulations for the training, validation, and test dataset as reported in Fig. 8(a), (b), and (c) respectively.
As can be observed in Fig. 8(a) and (b), for training and validation datasets, an R 2 value of 0.99 is obtained. For the test dataset, an R 2 value of 0.97 is obtained as shown in Fig. 8(c). Such high R 2 value for the unseen test dataset proves the excellent prediction efficacy of the ML-based approach presented here as the test dataset was completely kept hidden from the model during the training process. The forthcoming sub-section presents the SHAP analysis, to shed more light on the interpretations of the predicted results and to elucidate more on the relative importance of different mesoscale architectural features.

SHAP results for interpretation of the predicted responses
Models such as NN are considered black box models [72] due to their complexity and highly non-linear architecture, leading to low interpretability. Towards this, the NN model is integrated here with the SHAP technique [83] to demonstrate the interpretative ability of the presented approach. In SHAP, each feature is assigned with an importance value, and the influence of all the features on the predicted responses is evaluated. Fig. 9 demonstrates that the Poisson's ratio for cement mortar mesostructure is dominated by the volume fraction of voids (or porosity), followed by the aspect ratio, length of the major axis, and the number of voids present in the system. This indicates that the volume fraction of voids (or porosity) has the highest influence on the Poisson's ratio which can directly affect the aspect ratio and the length of the major axis. Such inferences can provide the strategies to enhance the auxetic response in cement mortar mesostructure. It is also interesting to observe that the number of voids along the Y-axis, which is also the loading direction, is relatively more sensitive than that along X-axis. ure, the color represents the feature value, and its corresponding xaxis (SHAP) value represents the contribution to the output Poisson's ratios.
Here, the blue color represents the low feature value, and the red color represents the high feature value. The SHAP value of zero represents the mean value of the Poisson's ratio from the dataset. From Fig. 10, it can be clearly observed that for the volume fraction of the voids, the positive SHAP value is associated with low volume fraction values. This indicates that the low volume fractions of voids increase the model's output (Poisson's ratios) to the right by up to 0.18 from the mean Poisson's ratio of 0.11 obtained for the entire database. Conversely, a high-volume fraction of the voids reduces the model's output by up to 0.4 from the mean value. A similar trend is observed for the aspect ratio, length of the major axis, as well as the number of voids along the Â and y-axis where low feature values increase the Poisson's ratio from the mean value and higher feature values, shift the Poisson's ratio to the left (Poisson's ratio reduces) from the mean by a varying amount. It is worth mentioning that some high feature of the number of voids along the Y-axis is associated with increasing the model's output. Such behavior can be considered as a mixed effect. However, due to a low number of such points, they can be regarded as outliers. Overall, it is observed that the input features with low feature value tend to draw the output Poisson's ratio closer to the positive mean value, whereas the higher feature values lower the Poisson's ratio from the positive mean value towards the negative domain inducing auxetic effect. Such interpretations allow an expert to ensure that the functions learned by the ML models are reasonable and physically sensible. While the SHAP algorithm depicts the expected influence of the chosen geometric parameters in inducing the auxetic behavior, it also sheds more light on the trend that the overall    volume fraction and the length of the major axis are the most influential parameters toward achieving auxetic response. The efficacy of SHAP implementation can be assessed by its ability to interpret such realistic trends. As such, the NN-based ML approach when coupled with the SHAP algorithm can not only confirm physicsbased phenomenon in the predicted responses but also can lead to trends that are not apparently visible. Fig. 11 shows the river flow plot for the Poisson's ratio using the NN model. While the violin plot focuses on the contribution of a given input feature toward different property values, the river flow plot focuses on the contribution of different input features towards a given property value.
The expected model value shown in the river flow plot corresponds to the model output when values of any of the input features are not available. In other words, when no information about the input features is available, the best estimate of the output value will be the mean of the dataset used for training the model. Adding the information about each of the input features one-by-one, followed by all possible combinations, provides insight into the individual and the collective roles played by each of the input features in governing the output property. To demonstrate this effect, each line in the SHAP river flow plot (see Fig. 11) corresponds to one data point. The line color represents the property value corresponding to the given data point, the red color is associated with high content values and blue with the low content values. The rise and fall of the line with respect to the expected value show how these respective features control the final property value. The features are arranged in increasing order of average SHAP value from left to right of the x-axis. For the volume fraction of the voids, it can be observed that as the volume fraction increases from lower to higher values, the Poisson's ratio drastically reduces from the expected mean value. In contrast to the volume fraction, the aspect ratio has a mixed effect where the trends are not uniform. For instance, for one particular combination, high content of the aspect ratio shows to increase the Poisson's ratio from the expected value whereas, for another combination, higher aspect ratio results in a decrease in the Poisson's ratio. This indicates that the trend with respect to the aspect ratio is governed by the other parameters in the model. A similar mixed trend is also observed for the length of the major axis and the number of voids along the Y-axis, whereas the number of voids along the X-axis is less significant towards the negative Poisson's ratio. This also indicates that the auxetic behavior is primarily controlled by the number of elliptical voids along the direction of the loading. Overall, the trends suggest that the low input features content favors the positive Poisson's ratio whereas the negative Poisson's ratio or auxetic response is highly influenced by the increasing volume fraction of the elliptical voids.

Conclusion
This paper presents an ML-based approach to predict the auxetic behavior of cellular cementitious composites enabled by judicious placement of elliptical voids in the mortar mesostructure. To address the unavailability of a large experimental dataset for such composites, this paper implements FE-based numerical simulations to generate a dataset of Poisson's ratios with 850 combinations by varying various structural features of the composite such as the volume fraction, aspect ratio, length of the major axis of the elliptical voids in the mortar mesostructure as well as the number of voids along x and y-axis. The dataset is leveraged to predict the auxetic behavior of the composites using a feed-forward multilayer perceptron-based NN approach. The optimal hyperparameters are judiciously selected where an optimal trade-off between predictive accuracy and the level of complexity is sought to avoid any under-or overfitting. The NN model with two hidden layers and 14 neurons yielded exceptional prediction efficacy. SHAP technique is also implemented for the interpretation of the responses predicted by the NN model, where the SHAP values for all the input parameters are calculated to shed light on their relative sensitivity towards the predicted Poisson's ratios. Based on SHAP values, it was found that the Poisson's ratio of the representative mesostructure of the mortar is most dominantly influenced by the volume fraction of the voids (or porosity) followed by the length of the major axis and the aspect ratio. Besides, the results also indicate that the number of elliptical voids along the direction of the loading shows a significantly higher influence on the auxetic behavior as compared to the number of voids perpendicular to the loading direction. Overall, the results in this paper demonstrate machine learning as a promising tool to predict the auxetic behavior of the cellular cementitious composites with judiciously placed elliptical voids, and as such, the ML-based predictive tool can enable materials designers and decision-makers, to make informed decisions regarding the selection of the design parameters to obtain the desired auxetic behavior in these cementitious composites for a wide variety of applications. Moreover, although this paper implements the ML-based predictive approach toward prediction of auxetic response in cementitious composites, the developed FE simulation-based dataset generation, and corresponding ML-based performance prediction approach are generic. As such, the approach can be extended toward a variety of structural metamaterials including a wide range of polymer composites by assigning appropriate material properties to the solid part to generate a large dataset and by optimizing the hyperparameters of the neural network accordingly.

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