Methodology for Optimizing Composite Design via Biological Pattern Generation Mechanisms

Mechanistic capabilities found in nature have influenced a variety of successful functional designs in engineering. However, the unique combinations of mechanical properties found in natural materials have not been readily adapted into synthetic materials, due to a disconnect between biological principles and engineering applications. Current biomimetic material approaches tend to involve mimicking nature's microstructure geometries or mimicking nature's adaptive design process through brute force element-by-element composite optimization techniques. While the adaptive approach promotes the generation of application-specific microstructure geometries, the element-by-element optimization techniques encompass a large design space that is directly related to the number of elements in the system. In contrast, a novel methodology is proposed in this paper that merges biological pattern generation mechanisms observed in the Gray-Scott model, with an evolutionary-inspired genetic algorithm to create adaptive bio-inspired composite geometries optimized for stiffness and toughness. The results reveal that this methodology significantly reduces the optimization parameter space from tens of thousands of parameters to only four or five. In addition, the resultant composite geometries improved upon the overall combination of stiffness and toughness by a factor of 13, when compared to published brute force element-by-element techniques.


Introduction
The design of many of mankind's technological advancements have been inspired by nature. Airplanes were inspired by the wings of eagles, the Japanese bullet train was inspired by the beak of the kingfisher, and Velcro was inspired by burrs [1]. These successes have come from mimicking the function of nature. The next technological breakthroughs will come from utilizing the design of natural materials to inspire the generation of novel synthetic structures and microstructures with enhanced properties and functionality for engineering applications. Natural materials often contain unique combinations of mechanical properties that are mutually exclusive in synthetic materials. For example, shell and bone are known to be strong and tough, enamel is known to be hard and wear resistant, and horns and antlers are known to be exceptionally damage tolerant [2][3][4]. Biomimetic research has revealed that the unique combination of mechanical properties can often be linked to the complex geometric microstructural architecture [5][6][7][8]. However, design concepts found in nature have not been readily adapted to synthetic materials, because it is still unclear how to apply bio-like principles to practical engineering applications [9]. One extensively researched bio-inspired approach involves mimicking structural design elements found in nature. Researchers have shown that nacre-inspired brick-and-mortar composite microstructures can increase a material's elastic modulus while also improving upon the toughness through crack deflection promoted by the geometric design [10][11][12]. Stomatopod-inspired helicoidal composite microstructures have been shown to increase energy absorbing capabilities due to a geometry that fosters longer cracking paths [12][13][14]. Conch shell-inspired cross-lamellar composite microstructures have been shown to enhance impact resistance and stiffness [12,15]. Finally, osteon-inspired tubular composite microstructures have been shown to increase the strength and toughness of the material through a wide array of toughening mechanisms [12,16]. The primary limitation of this approach is that a majority of the bio-like geometries are generated using rigid building blocks that are repeated throughout the structure. Typically, these building blocks contain overly simplified versions of complex microstructural architectures found in nature. While the simplified geometries produce enhanced material properties, they lack the ability to be tailored to specific engineering applications. Therefore, while the approach may capture some of the basic geometries found in nature's composites, it lacks the critical adaptive component and geometric complexity of the design processes found in nature.
Rather than directly mimicking nature, another design approach involves generating optimized composite architectures via a quasiconvex variational approach. The theory is based on the idea that a composite is optimal if the stress field is evenly distributed across the material domain, thereby minimizing the total strain energy [17][18][19]. While this definition of an optimal composite is still widely accepted, the methodology of generating the microstructure designs continues to evolve. For example, a more recent bio-inspired approach attempts to mimic the adaptive processes found in nature by designing composite microstructures, on an element-by-element basis, for a particular loading configuration. This approach permits the generation of customizable composite architectures that can be tailored to meet both the local and global demands placed on the material. Current research in this area of biomimetics has primarily focused on optimizing the distribution of stiff and compliant material constituents across a gridded domain [20][21][22][23][24].
The work of Gu et al. and Abueidda et al. present case studies for optimizing the stiffness and toughness of a pre-cracked composite plate specimen [20,21]. In the work of Gu et al., a modified greedy algorithm was used to evaluate the material assigned to each element in a 20 × 20 grid and determine if it should be stiff or compliant [20]. In the work of Abueidda et al., a genetic algorithm was used to optimize a 16 × 16 grid by encoding the material assigned to each element into a set of genes describing the individual chromosome [21]. In order to reduce the overall size of the problem, both studies implemented a plane of symmetry along the length of the crack; additionally, Gu et al. placed a 20% volume fraction restriction on the compliant material to further reduce the size of the problem [20,21]. These studies revealed that their adaptive approaches were capable of evolving a random solution into a composite microstructure optimized for stiffness and toughness. However, the primary limitation of this approach is that the optimization process must consider every element in the domain as a tunable parameter. For a symmetric 20 × 20 grid, 200 parameters have to be evaluated each time to produce a new composite geometry. If the symmetry conditions are removed, and/or the grid size becomes more refined, and/or the number of dimensions increase, the number of optimization parameters also increase.
This work aims to introduce a new adaptive approach to the design of bio-like structures by implementing a bio-inspired pattern generation algorithm to generate the composite microstructure geometries, rather than designing them on an element-by-element basis as previous researchers have done [20,21]. The proposed algorithm employs biological mechanisms that have been used to study morphogenesis and pattern development in nature, in conjunction with an evolutionary algorithm. It is capable of producing some of the most abundant patterns found in nature such as spots, stripes, and hexagons. The category of pattern produced and the refinement of the pattern within a category are controlled using four tunable parameter settings embedded within the algorithm. As a result, regardless of how finely resolved the domain becomes, the microstructure geometries can be optimized by varying only four parameters. The proposed methodology is shown to produce significantly improved application-specific microstructures relative to the previous element-by-element attempts.

Methodology
The composite microstructure geometries, designed with the bioinspired pattern generation algorithm, were analyzed using the same material constituents and loading configuration employed in the work of Gu et al. and Abueidda et al. [20,21]. This permitted the direct comparison to a set of known results. The loading configuration was comprised of a pre-cracked plate specimen under uniaxial tension as illustrated in Fig. 1. The initial crack length was defined as a percentage of the total edge length. The crack was centered along the specimen's edge and displacement boundary conditions were applied in the x-direction perpendicular to the length of the crack. The materials used in the design of the composite microstructures were defined as stiff or compliant. The objective of this study is to examine how well the bio-inspired composite microstructures can adapt to the applied loading configuration when optimized for stiffness and toughness.

A. Gray-Scott model
The bio-like composite microstructure geometries were generated with the Gray-Scott model. It was selected due to its unique ability to produce a wide array of patterns found in nature using biological pattern generation mechanisms. Initially, the model was created by Gray a w Fig. 1. A pre-cracked square plate specimen with a crack length of a = 0.2w or a = 0.25w placed under a uniaxial tension test with displacement boundary conditions in the x-direction. and Scott to study the autocatalytic chemical reactions of an isothermal system in a continuously fed, well-stirred, open-tank reactor [25][26][27]. However, Pearson extended the model into two-dimensions and showed that a variety of chaotic and Turing patterns evolve with time and space in response to the perturbation of a steady-state system [28]. The reaction-diffusion equations he developed are given by where U and V represent the concentration values of chemicals U and V respectively, D u and D v are the diffusion coefficients of the chemicals, F is the feed rate, and k is the kill rate [28]. Within the set of equations, the feed rate and kill rate control the type of pattern that is produced while the diffusion coefficients control the length scale. The generation of patterns via reaction-diffusion kinetics had been previously studied in Turing's work on morphogenesis (the study of pattern formation in nature) in which patterns such as spots, stripes, hexagons, and spirals became termed "Turing patterns" [29]. As a result, a majority of previous studies involving the Gray-Scott model have primarily focused on pattern characterization and development [30][31][32][33][34].
To date, the patterns produced by the Gray-Scott model have not been analyzed as load bearing composite materials. The Gray-Scott model utilizes diffusion mechanisms between two chemical species as the driving force in pattern formation. Therefore, the patterns are computationally represented as two continuous variables, given by the concentration values of chemical U and chemical V, across a gridded domain as illustrated in Fig. 2a. However, in order to convert the patterns into bi-material composite microstructures, the continuous domain must be discretized into a binary representation. A previous study measuring the complexity of the Gray-Scott model binarized the gridded domain by assigning elements a value of 1 if the concentration value of chemical V exceeded 0.3 and a value of 0 otherwise [32]. While the approach is simple, it is incorrect to assume that all patterns produced by the Gray-Scott model have the same constant cutoff value. In fact, by making this assumption, there is a possibility that some of the patterns that show up in the continuous domain become homogeneous in the discrete domain.
A bifurcation diagram produced by Mazin et al. reveals a saddlenode bifurcation curve that determines the stability of the chemical system [35]. The curve is given by the equation the concentration value of chemical V and F is the feed rate; above this curve the system may be stable whereas below the curve the system is always unstable unless the concentration value of chemical V is zero [35]. Using the saddle-node bifurcation curve, the gridded domain can be discretized into a binary representation by assigning elements a value of 0 if the concentration value of chemical V is less than the square root of the feed rate and a value of 1 otherwise, as described in Eq. (3) [36].
Since the patterns are generated using a range of feed rates from 0 to 0.25, the cutoff values for discretizing the domain vary from 0 to 0.5, depending on the type of pattern that was generated. As a result, this discretization technique ensures that the stable portion of the patterns are always captured regardless of the type of pattern produced. Fig. 2 illustrates this technique by converting a spotted pattern produced by the Gray-Scott model in the continuous domain to the discrete domain. Fig. 2a displays the continuous domain as a gray-scale image. Fig. 2b reveals the 3D surface plot of the concentration values of chemical V across the gridded domain. The black cutoff plane was generated using the pattern's respective feed rate and the equation V ¼ ffiffiffi F p . Using the information provided in Fig. 2b, the patterns were converted to the discrete domain by assigning elements above the plane a value of 1 and elements below the plane a value of 0. Once the patterns were binarized, they were converted to composite microstructures by assigning the 0's to the stiff material constituent and the 1's to the compliant material constituent, as illustrated in Fig. 2c.
Due to the reaction-diffusion kinetics that evolve the patterns in time and space, careful attention must be paid to the timestep, stopping criteria, and seeding location. If the time step is too large, the concentration values in the grid can become negative or undefined numerical values. If the time step is too small, the speed of the simulation is degraded. According to Mazin et al., numerical stability is generally achieved with a maximum stability index of 0.4 [35]. Therefore, in order to ensure numerical stability, a stability index of 0.2 was enforced. Due to the variability in diffusion coefficients, each pattern required a unique time step for maximum simulation speed as defined by Eq. (4). Within Eq. (4), Δt is the time step, 0.2 is the stability index, Δx is the spatial resolution, and D u is the diffusion coefficient of chemical U. However, higher ratios of diffusion coefficients can require reduced stability indexes. Therefore, the initial time step described in Eq. (4) was decreased by a factor of 1.125 each time a negative or undefined numerical value occurred.
The stopping criteria are a critical component of the Gray-Scott model as some classes of patterns (such as chaotic patterns) are unstable and can evolve indefinitely. In order to prevent this from occurring, the patterns were allowed to evolve until they met one of the following stopping criteria: 1. The predefined volume fraction for the stiff material was met. 2. The pattern became homogeneous as indicated by a volume fraction of one or zero. 3. The propagating pattern reached all four boundaries of the grid. At this point in time, the pattern was allowed to propagate for an additional 5000 iterations before ceasing. 4. The change in volume fraction was less than 1% over a period of 5000 iterations. 5. The number of iterations exceeded 100,000. These multiple stopping criteria were implemented to ensure the evolution of time sensitive solutions.
The final consideration for the Gray-Scott model is the seeding location, or the area that the pattern propagates from. Two different seeding locations were tested to explore the impact on the final microstructure geometry and resultant composite material properties. Due to the presence of an edge crack in the specimen, the first seeding location was defined as a rectangular area surrounding the crack, with a width of 20 elements and a height of 46 elements, as illustrated by the blue striped region in Fig. 3a. The second seeding location was located along the specimen's cracked edge, with a height of 10 elements and a 10 element spacing from each outer border, as illustrated by the blue stiped region in Fig. 3c. The grid was initially placed into a steady-state environment with U = 1 and V = 0. The seeding location was then perturbed to U = 0.5 and V = 0.25, thereby initiating pattern propagation of the compliant material [28]. Fig. 3b and d reveal how various patterns propagate from their respective seeding locations.
The patterns were generated with a spatial resolution of Δx ≅ 0.01 on a 144 × 144 square mesh grid. The kill rate was allowed to vary from 0 to 0.1 and the feed rate was allowed to vary from 0 to 0.25, thereby encompassing the full pattern generation region. Forward Euler integration, Neumann boundary conditions, and the five-point finite-difference Laplacian approximation were implemented. Deviating slightly from the work of Pearson, the diffusion coefficients were also allowed to be varied rather than held at a constant D u = 2 * 10 −5 and D v = 1 * 10 −5 [28]. Mazin et al. mathematically proved that increasing the ratio of diffusion coefficients, σ = D u /D v , also increased the area where Turing patterns could be generated [35]. Therefore, to maximize the pattern generation region, the diffusion coefficients were allowed to be varied as long as the diffusion coefficient of chemical U was always greater than the diffusion coefficient of chemical V. In our previous work, we also discovered that varying the diffusion coefficients permits the generation of a new class of pattern defined by ν in Fig. 4 [36]. As far as the authors know, this new class of pattern has not been recorded in any of the literature on the Gray-Scott model. Patterns α through μ in Fig. 4 were generated from Pearson's work, while patterns ν were generated by varying the diffusion coefficients [28]. The set of 14 patterns and their corresponding parameter settings, as illustrated in Fig. 4, were implemented as the initial population in the genetic algorithm described in Section 1.4. set to 25% of the specimen edge length when comparing to the work of Abueidda et al. and 20% of the specimen edge length when comparing to the work of Gu et al. [20,21]. In order to match the resolution of the Gray-Scott model, a 144 × 144 square mesh grid, with plane stress 4-node bilinear reduced integration elements (CPS4R) and two degrees of freedom, was implemented. Since the patterns produced by the Gray-Scott model were converted into microstructures with two definable regions, each region could be assigned its own set of material properties. The black region was characterized as a stiff material with a Young's modulus of E s = 1GPa, a Poisson's ratio of ν s ¼ 1 3 , and a strain to failure of ε fs = 0.1. The white region was characterized as a compliant material with a Young's modulus of E c = 0.01GPa, a Poisson's ratio of ν c ¼ 1 3 , and a strain to failure of ε fc = 1. The stiff and compliant materials were assumed to be isotropic linear elastic. The applied boundary conditions consisted of a small uniaxial displacement in the x-direction, perpendicular to the length of the crack, as shown in Fig. 1.

Finite element analysis
Upon completion of the simulation, the maximum principal stress and strain for each element were output from the model. Based on this information, the stiffness and toughness of the overall composite microstructure could be calculated. The stiffness was defined as the slope of the stress-strain curve or the Young's modulus of the composite structure. It was calculated using the average stress and strain across all of the elements as defined by Eq. (5); where E is the Young's modulus of the overall composite microstructure, σ AVG is the average stress in the loading direction across all of the elements, and ε AVG is the average strain in the loading direction across all of the elements. The toughness was defined as the energy required to initiate crack growth. As a result, the composite was considered failed once the crack began to propagate, or the strain at the crack tip reached the failure strain of the material. Due to the linear elastic behavior of the materials and the potential for non-symmetric microstructure geometries, a failure index for each of the two elements at the crack tip, as described in Eq. (6), was calculated to determine which element would fail first. The maximum failure index was then used to calculate the toughness of the overall composite microstructure as defined by Eq. (7). Within Eqs. (6) and (7), i = L is the element left of the crack tip, i = R is the element right of the crack tip, F i is the failure index, ε tip(i) is the current strain in the element, ε f(i) is the strain to failure associated with the material assigned to that element, and T tip is the toughness of the composite structure as defined by the crack tip elements.
In contrast to the work of Gu et al. and Abueidda et al., this paper also investigates a broader definition of toughness that is based on the failure index of all of the elements in the composite rather than focusing solely on the elements at the crack tip [20,21]. Gu et al. revealed that surrounding the crack with the compliant material can reduce the strain at the crack tip, thereby increasing the crack tip definition of toughness [20]. While this definition of toughness is mechanically correct, it does not consider the possibility of other elements in the composite failing before the crack tip. Therefore, this paper introduces a broader definition of toughness that deems the composite failed when the first element, regardless of its location, reaches the failure strain of the material. In order to calculate this, a failure index was generated for every element in the composite and the maximum failure index was used to calculate the toughness according to Eqs. (8) and (9); Fig. 4. Set of 14 patterns, generated with the Gray-Scott model, and their corresponding parameter settings (feed rate (F), kill rate (k), diffusion coefficient of chemical U (D u ), and diffusion rate of chemical V (D v )) that were implemented as the initial population in the genetic algorithm. Patterns α through μ were generated from Pearson's work and converted into two component microstructures (black and white) [28]. The patterns labeled ν appear when the diffusion coefficients are allowed to be varied.
where j corresponds to the element number and T comp corresponds to the toughness of the composite structure as defined by the element with the maximum failure index.

Genetic algorithm
A genetic algorithm was implemented to optimize the bio-like composite microstructure geometries, generated with the Gray-Scott model, for stiffness and toughness. Three different case studies were analyzed based on the current research being conducted in the field of adaptive biomimetic composite materials. A subscript of one indicates that the equation or material property values were adopted from the work of Gu et al. [20]. Similarly, a subscript of two indicates that the equation or material property values were adopted from the work of Abueidda et al. [21]. Each case study utilized a different fitness function to quantify the performance of the various microstructure geometries, in terms of the composite's overall stiffness and toughness. The fitness function is essential for optimization in a genetic algorithm, as it helps guide the search through the parameter space. As a result, the genetic algorithm's fitness function employed in each case study is detailed in the following paragraphs and summarized in Table 1.
The first case study was based on the work Gu et al. who optimized the stiffness and toughness of the composite materials on an elementby-element basis using a modified greedy algorithm [20]. In order to stay consistent with their work, the pre-cracked plate specimen had a crack length equal to 20% of the edge length, the volume fraction of the stiff material was restricted to 80%, and the objective function for the greedy algorithm was implemented as the fitness function in the genetic algorithm. In the work of Gu et al. the objective function was defined according to Eq. (10), where E is the Young's modulus, T is the toughness, and E max(1) and T max(1) are the normalization parameters that were pre-selected from the initial population [20]. For comparative purposes, the optimal composite geometry (for maximum stiffness and toughness) obtained from the work of Gu et al. was run using the FEA model developed in this study to acquire the resultant elastic modulus and toughness which were implemented as the normalization parameters [20]. Therefore, E ð Þ ¼ 1 indicates that the microstructures produced by the Gray-Scott model had the same elastic modulus or toughness values as the optimal solution produced in the work of Gu et al. [20]. Finally, in the event that a microstructure generated by the Gray-Scott model was unable to meet the required 80% volume fraction of the stiff material, the fitness function was modified to penalize the result, as illustrated in Eq. (10).
The second case study was based the work of Abueidda et al. who optimized the stiffness and toughness of the composite materials on an element-by-element basis using a genetic algorithm [21]. In order to stay consistent with their work, the pre-cracked plate specimen had a crack length equal to 25% of the edge length and no pre-existing restrictions were placed on the volume fraction of the materials. In the work of Abueidda et al., the fitness function was defined according to Eq. (12) where E is the Young's modulus, T is the toughness, and E max(2) and T max (2) are the normalization parameters that were acquired from single parameter optimization of the stiffness and toughness [21]. It's important to note that Abueidda et al. utilized a neural network in conjunction with the genetic algorithm to develop their optimal solutions [21]. However, the neural network introduces error which could affect the final results. Therefore, the multi-parameter optimal solution (for maximum stiffness and toughness) was selected for comparative purposes and run using the FEA model developed in this study, to acquire the normalized elastic modulus and toughness parameters [21]. Furthermore, ð Þ ¼ 1 indicates that the microstructures produced by the Gray-Scott model had the same elastic modulus or toughness values as the optimal solution produced in the work of Abueidda et al. [21]. Finally, while no restrictions were placed on the volume fraction of the stiff material, the Gray-Scott model required a pre-defined volume fraction for its primary stopping criteria. As a result, the genetic algorithm was allowed to optimize the volume fraction of the stiff material by incorporating the parameter as gene in the individual chromosomes.
The third case study explored the implementation of the broader definition of toughness defined in Eq. (9). The composite specimen was defined with a 25% crack length and no restrictions were placed on the final volume fraction of the stiff material. The genetic algorithm utilized the fitness function defined in Eq. (10). However, the normalization parameters E max and T max were defined by the geometry with the maximum stiffness and toughness from the initial population of microstructures generated with the Gray-Scott model. Table 1 summarizes the three different case studies presented in this work.
While the fitness function depicts the performance of the overall composite microstructure geometry, it must also guide the genetic algorithm during the optimization process. The genetic algorithm mimics various adaptation mechanisms implemented by organisms in natural systems to produce superior solutions [37]. As a result, the algorithm is comprised of a population that changes from generation to generation. This population is made up of a variety of individuals referred to as chromosomes. Each individual chromosome is made up of a set of genes that describe the characteristics of the individual. In this study, the genes contained information about how to generate a pattern using the Gray-Scott model. In the first case study, each individual chromosome contained four genes comprised of the Gray-Scott model's four parameter settings: feed rate (F), kill rate (k), diffusion coefficient of chemical U (D U ), and the diffusion coefficient of chemical V (D V ). In the second and third case studies, each individual chromosome contained five genes comprised of the Gray-Scott model's parameter settings and a pre-defined volume fraction of the stiff material (VF). Since there was no pre-existing restriction on this value, as in the first case study, the additional gene was implemented to optimize the Table 1 Details on the three different case studies that were studied based on the current research being conducted in the field of adaptive biomimetic composite materials. Shown below is the case study, citing work, crack length in terms of the length of the specimen's edge, restrictions placed on the volume fraction of the stiff material, the definition of toughness that was implemented, the fitness function used and its respective equation number, and the parameters that were optimized in the genes of the genetic algorithm.  Table 2.
Once the chromosomes are defined, they must undergo various adaptation procedures aimed at improving the overall fitness of the population from one generation to the next. The genetic algorithm implemented in this study utilized the traditional genetic operators of selection, crossover, mutation, and elitism. The selection operator is used to select the parent chromosomes that will breed the next generation. It plays a key role in how effectively the genetic algorithm converges to a global optimum by maintaining a certain degree of diversity within the population [38]. The crossover and mutation operators are used to balance the amount of exploration and exploitation across the parameter space. Exploration refers to a global search of the parameter space into regions where the solution is not well defined. Exploitation refers to a local search of the parameter space where the goal is to improve upon pre-existing well-known solutions. In terms of a genetic algorithm, exploitation is generally achieved through the crossover operator, while exploration is generally achieved through the mutation operator [39]. Finally, the elitism operator is used to eradicate any degradation of the solution space due to the stochastic process. Elitism works by passing the fittest chromosomes from one generation to the next [40].
During the selection stage, the roulette wheel technique was implemented to give each individual chromosome a certain probability of being selected. Chromosomes with higher fitness values were favored by a greater probability of being selected than chromosomes with a lower fitness value. Since the chromosomes are made up of a relatively small number of genes, the diversity of the population was quickly degraded. In order to prevent this from occurring, four sets of parents were selected to breed 12 chromosomes each and two unique chromosomes with the highest fitness parameters from the current generation were always carried into the following generation through the elitism operator. As a result, each generation had a total population size of 50 chromosomes. In addition, the diversity of the chromosomes in the first generation also played a key role in the diversity of the following generations. Therefore, the first population was generated by applying the selection, crossover, mutation, and elitism operators to the set of patterns described in Figure 4 that represent the various classes of patterns that the Gray-Scott model is known to produce; thereby ensuring that the first population was comprised of a diverse selection of chromosomes.
Furthermore, during the selection process, the diversity was degraded due the selection of homogeneous patterns that could shift the population into a region of the parameter space where no patterns are produced. In order to avoid this scenario, the microstructures were evaluated using two rapid analysis techniques termed morphology movement and connectivity. The morphology movement (MM) utilizes the binarized matrix (BM) to subtract sequential layers from one another in a single direction, and sum the number of differing elements from one layer to the next; as illustrated in Eq. (13) where n is the total number of elements in the desired direction.
This technique provides a quantitative representation of the amount of change that occurs throughout the composite's microstructure. If a pattern is homogeneous, the morphology movement in both the x and y directions will be relatively small, as there are no complex movements exhibited by the microstructure. In addition, a homogeneous pattern can be identified as a single structure that is fully connected. As a result, the homogenous patterns were identified by analyzing the combination of the microstructure's morphology movement and connectivity. Once identified as a homogeneous structure, it was assigned a fitness of zero to ensure that it couldn't be selected and undermine the diversity of the population.
During the crossover and mutation stages, a crossover probability of 90% and a mutation probability of 20% were implemented. The single point crossover technique was used to swap the genes between the parent chromosomes. Due to the continuous nature of the genes, mutation was performed using a lognormal distribution in order to ensure the generation of solely positive values. The distribution was centered around the value of the original gene so that the mutated solution would have a high probability of being close to the original solution. If the mutated solution exceeded the bounds placed on the gene, remutation was performed until compliance was met. Finally, the genetic algorithm was allowed to run for 120 generations or until there was no significant change in the maximum fitness value, at which point the optimization process was terminated. A summary of the genetic algorithm's parameter settings can be found in Table 2.

Results
For each case study, the composite microstructure geometries produced by the Gray-Scott model were optimized for stiffness and toughness. In addition, the Gray-Scott model was seeded from two different locations (an area surrounding the crack and an area along the cracked edge) in order to examine the effects on the final microstructure geometry and the overall fitness of the composite material. The top three structures with the highest fitness values for each case study and seeding location are presented. These results were obtained from running the genetic algorithm multiple times to ensure a thorough search of the parameter space. Multiple simulations act to increase the confidence that an optimal solution was obtained; since the genetic algorithm does not ensure convergence of the global optimum.
The top microstructure geometries with the highest fitness values from the first case study are presented in Fig. 5a-5b. Each composite geometry had a stiff material volume fraction of 80%. Fig. 5a displays the top three microstructure geometries generated from the crack seeding location with fitness values of 13.3, 10.1, and 3.8 from top to bottom respectively. Fig. 5b displays the top three microstructure geometries generated from the edge seeding location with fitness values of 9.0, 7.6, and 5.9 from top to bottom respectively. Ultimately, the crack seeding location produced the geometry with the highest overall fitness value. However, regardless of the seeding location, all of the microstructure geometries had a fitness value greater than one. These values indicate that the methodology established in this paper was capable of producing geometries that exceeded the optimal solution presented in the work of Gu et al. [20]. Fig. 5c displays the normalized stiffness and toughness values for each microstructure geometry. The plot reveals that on average the crack seeding location produced composite geometries that were better suited for maximizing toughness, while the edge seeding location produced composite geometries that were better suited for maximizing stiffness. Finally, Fig. 5d reveals how the best solution from each generation of the genetic algorithm evolved throughout the optimization process for each seeding location. The top microstructure geometry for the crack seeding location reached its maximum fitness value (13.3) after 108 generations of the genetic algorithm, while top microstructure geometry for the edge seeding location reached its maximum fitness value (9.0) after 111 generations of the genetic algorithm. From this plot it is clear that the genetic algorithm was able to successfully tune the Gray-Scott model's parameter settings towards improved microstructure geometries. The horizontal lines in the plot occur due to the implementation of the elitism operator which prevented any degradation of the maximum fitness value. The top microstructure geometries with the highest fitness values from the second case study are presented in Fig. 6a-6b. Only one geometry from each seeding location is represented, because they were the only geometries that the model converged to. Fig. 6d illustrates how the best microstructure geometries from each generation evolved throughout the optimization process for multiple runs of the genetic algorithm. This plot illustrates that the genetic algorithm suffered from premature convergence. The microstructure geometry for the crack seeding location reached its maximum fitness value (1.8) after an average of 37 generations of the genetic algorithm, while the microstructure geometry for the edge seeding location reached its maximum fitness value (2.3) after an average of 25 generations of the genetic algorithm. Fig. 6c reveals that the premature convergence occurred due to the fitness function informing the genetic algorithm to solely focus on optimizing the stiffness. These results highlight the importance of a well-balanced fitness function that reduces the fitness if one material property is significantly less than the other. As a result, when optimizing the Gray-Scott model's parameter space, careful attention must be paid to the general form of the fitness function.
The top microstructure geometries with the highest fitness values from the third case study are presented in Fig. 7a-7b. Fig. 7a displays the top three microstructure geometries generated from the crack seeding location with fitness values (and stiff material volume fractions) of 4.0 (86%), 3.95 (87%), and 3.88 (90%) from top to bottom respectively. Fig. 7b displays the top three microstructure geometries generated from the edge seeding location with fitness values (and stiff material volume fractions) of 14.6 (79%), 14.5 (81%), and 14.4 (86%) from top to bottom respectively. In order to compare the fitness values from the first and third case studies, the fitness value of the top microstructure geometry from the first case study (Fig. 5a with a fitness value of 13.3) was recalculated using the fitness function and broader definition of toughness implemented in the third case study. This resulted in a newly calculated fitness value of 2.2. The broader definition drove down the overall value of toughness by approximately three orders of magnitude, due to failure occurring elsewhere in the composite other than at the tip of the crack, as illustrated by the failure index map in Fig. 8a. Therefore, the optimal composite geometry for the crack tip definition of toughness was the worst composite geometry for the broader definition of toughness. Ultimately, the broader definition of toughness revealed that, for this particular loading configuration, complex microstructure designs are not necessary, as illustrated by the top geometries in Fig. 7b. It also revealed that the Gray-Scott model was capable of adapting its solutions, from simple to complex topologies, based on the requirements of the system. In addition, the third case study permitted the optimization of the volume fraction of the stiff material. In general, the volume fraction was greater than the 80% restriction placed on the geometries in the work of Gu et al. [20]. Therefore, in order to obtain the most optimized geometries possible, the volume fraction must be an optimizable parameter setting. Fig. 7c displays the normalized stiffness and toughness values for each microstructure geometry. This plot illustrates that the edge seeding location promoted the generation of composite geometries that were better suited for maximizing both stiffness and toughness. Finally, Fig. 7d reveals how the best solution from each generation of the genetic algorithm evolved throughout the optimization process for each seeding location. The top microstructure geometry for the crack seeding location reached its maximum fitness value (4.0) after 95 generations of the genetic algorithm, while top microstructure geometry for the edge seeding location reached its maximum fitness value (14.6) after 50 generations of the genetic algorithm. From this plot, it is clear that the genetic algorithm was able to successfully tune the Gray-Scott model's parameter settings and stiff material volume fraction towards microstructure geometries with increased fitness values.

Discussion
The authors note that the grid used in this study was significantly more refined than the grid implemented in the work of Gu et al. and Abueidda et al. [20,21]. As a result, the ability to generate more intricate microstructures could have aided in the development of geometries with increased stiffness and toughness. However, the improvement in the material properties was also likely due to the substantial reduction in the overall size of the design space that was being optimized. For example, if Gu et al. and Abueidda et al. were to optimize the 144 × 144 grid with no symmetry conditions enforced, their element-by-element techniques would require the optimization of 20,736 parameters [20,21]. On the contrary, by implementing the methodology presented in this paper, the microstructure geometries generated with the Gray-Scott model were optimized by tuning only four or five parameters. In addition, if the grid were to become even more refined and/or the number of dimensions were to increase, the number of parameters required in the element-by-element method would grow proportionally, while the number of parameters required with the methodology presented in this paper would remain at four or five. Therefore, by using biological mechanisms to generate the composite microstructures, rather than designing them on an element-by-element basis, leads to a significantly smaller design space that could aid in the development of geometries with superior mechanical properties.
In terms of evaluating the mechanical properties, two definitions of toughness were implemented throughout the various case studies. The first and second case studies were performed by defining toughness as the energy required to initiate crack growth. Therefore, the composite was not considered failed until the crack began to propagate. While this definition of toughness is satisfactory for a homogeneous material, it is not broad enough for composite materials. For example, by surrounding the crack with the compliant material, it is likely that failure will occur elsewhere in the composite before there is failure at the  crack tip. In order to illustrate this, the failure index of each element in the composite was mapped for the top crack and edge seeding location geometries from the first and third case studies, as displayed in Fig. 8. Fig. 8a-8b presents the failure index maps for the top geometries in the first case study. From these maps, it is evident that the crack tip definition of toughness was properly exploited, because the failure index at the crack tip was minimized. However, these maps also reveal that failure would have occurred in other elements first before there would have been failure at the crack tip; as indicated by the dark red elements which symbolize 100% failure. Fig. 8c-8d presents the failure index maps for the top geometries in the third case study where the broader defintion of toughness was implemented. Visually, the primary difference between the failure index maps is that the crack tip defintion of toughness had concentrated areas of high failure indicies, while the broader defintion of toughness had more dispersed areas of high failure indicies. As a result, the dispersion of failure decreased the maximum failure index (prior to normailization) for the crack seeding location from 52% for the crack tip definition of toughness (Fig. 8a) to 48% for the broader definition of toughness (Fig. 8c). Similairly, it decreased the maximum failure index (prior to normailization) for the edge seeding location from 67% for the crack tip definition of toughness (Fig. 8b) to 33% for the broader definition of toughness (Fig. 8d). Therefore, the broader definition of toughness promoted the generation of microstructure geometries that were ultimately further from failure than the geometries produced using the crack tip definition of toughness.

Future work
The novel adaptive methodology presented in this work was capable of evolving bio-like geometries, created using natural pattern generation mechanisms, into composite microstructures optimized for stiffness and toughness. This capability can readily be extended to encompass much more complex structures. In this study, the microstructure geometries were generated and optimized at a single length scale. However, it is well known that nature often utilizes multiscale hierarchical architectures to enhance the mechanical properties [8]. One of the major challenges facing the field of biomimetics is the lack of a methodology that can design and control multiscale hierarchical microstructures to improve mechanical properties for use in specific applications [9]. As a result, Fig. 9 presents an extension of this work that could potentially disrupt the field of biomimetics through controllable hierarchical microstructure designs.
The hierarchical microstructures illustrated in Fig. 9 were generated using the Gray-Scott model. The length scale was controlled through variable diffusion coefficients, while the geometry was controlled by adjusting the model's feed and kill rates. First, the macrostructure geometries were allowed to evolve until they met a desired stopping criterion, such as those presented in this paper. Next, the boundaries between the stiff and compliant material constituents were defined using Dirichlet boundary conditions; thereby constraining the pattern generation region, for the next level of hierarchical structure, to the region inside the macrostructure. Finally, in order to generate the next   level of hierarchical structures, the grid was refined, and the diffusion coefficients were decreased by an order of magnitude of 10 or more. Therefore, optimization of the hierarchical microstructure geometries could be performed by varying as few as eight parameter settings. For future work, we plan to extend the methodology presented in this study to encompass the generation of adaptive multiscale biolike composite architectures.

Conclusion
In summary, this paper presented a new paradigm in the field of adaptive biomimetic composite materials. This work introduced a novel methodology that utilized biological pattern generation mechanisms to create the bio-like composite microstructure geometries rather than designing them with the brute force element-by-element techniques. Thus, the design space during the optimization process was significantly reduced from 20,736 parameters for a 144 × 144 grid to only four or five parameters. As a result, this methodology permitted the generation and optimization of more intricate microstructures that ultimately led to an improved combination of stiffness and toughness by a factor of 13 over previous approaches. In addition, the proposed methodology presents a unique way of merging the biological processes observed in nature into the design of adaptive synthetic composite structures and microstructures. In nature it is not the full blueprint that gets stored in an organism's genes, rather it is a set of rules that can be modified based on the requirements of the surrounding environment [41]. In a similar matter, the methodology mimics this biological principle by optimizing a pattern generation process to produce the composite microstructure geometries.