Optimization of a breeding program design using evolutionary algorithms

Abstract

The present invention is in the field of breeding program design and concerns computer-implemented methods for the optimization of a breeding program design, comprising the steps of: (1) definition of the multivariate optimization problem and at least one termination criterion; (2) initialization of the first set of parameterizations; (3) evaluation of the first set of parameterizations based on the optimization problem defined in step (1); (4) selection of a number of parameterizations from the first, second or further set of parameterizations based on the evaluation made in step (3) or (6), wherein said number of selected parameterizations is lower than the total number of parameterizations in the set of parameterizations; (5) generation of a second or further set of parameterizations based on the selected parameterizations of step (4) comprising (5a) creating new parameterizations from the selected set of parameterizations obtained in step (4) by combination of selected parameterizations; and / or (5b) creating new parameterizations from the selected set of parameterizations obtained in step (4) by modification of selected parameterizations; (6) evaluation of the second or further set of parameterizations generated in step (5) based on the optimization problem defined in step (1) and derivation of optima, wherein unless at least one termination criterion is met steps (4) to (6) are repeated; and (7) assessment of the optima derived in step (6). The invention further relates to a data processing system comprising means for carrying out said methods, a computer program comprising instructions to if executed by a computer allow carrying out the methods, and a data carrier having stored thereon said computer program.

Description

[0001] OPTIMIZATION OF A BREEDING PROGRAM DESIGN USING EVOLUTIONARY ALGORITHMS

[0002] Field of the invention

[0003] The present invention is in the field of breeding program design. It is based on the development of a new optimization framework that makes use of stochastic simulation to evaluate the potential of specific breeding programs and integrates the use of evolutionary algorithms to iteratively optimize breeding programs.

[0004] Background of the invention

[0005] With the rise of genomics, advancements in biotechnology, statistical modeling, and shifts in market demands, breeding programs have undergone a substantial transformation in recent decades. The strategic combination of these advancements empowers breeders to intricately refine and adapt their breeding strategies. As a result, modern breeding programs are a complex resource allocation problem to manage both genetic gain and long-term sustainability

[0006] Breeding programs require substantial investments in both resources and time. Moreover, the consequences of breeding decisions may only be evident after several years. In addition to operating costs and complex logistics associated with it, breeders must consider a variety of design parameters related to specific breeding objectives, taking into account the potential risks and uncertainties associated with each decision and weighing the costs and benefits of each possible resource allocation

[0007] This is further complicated by the fact that breeding actions or changes to breeding program parameters are highly interdependent, and a change in one step of the breeding program will impact a multitude of key characteristics (e.g., genetic gain, genetic diversity, cost) of the breeding program.

[0008] Therefore, breeders seek to understand and assess the uncertainties that are inherent in predicting the outcomes of their breeding programs, given the long-term nature and complexity of the breeding process

[0009] A strategy that recently gained in popularity is to use stochastic simulation to assess potential changes or modifications in existing breeding program design before practically implementing them or using them as a tool for designing new breeding programs.

[0010] However, the optimization of a breeding program design using stochastic simulation is complicated by the fact that the output of a simulation is only the realization of a stochastic process. Thus, multiple replicates are necessary to reliably estimate the expected outcomes of a breeding scheme.

[0011] Since breeding programs often involve numerous parameters, it is not feasible to simulate all possible breeding designs many times, as simulating a real-world breeding scheme is computationally expensive and time consuming. Therefore, analysis of breeding program designs using stochastic simulations is usually limited to a couple of potentially interesting scenarios and research studies focusing on very specific aspects of breeding design. Recently, a framework for optimizing breeding program designs to address and generalize different aspects of the breeding program more effectively has been reported by Hassanpour et al. (Hassanpour et al. (2023) Optimization of breeding program design through stochastic simulation with kernel regression. G3 Genes|Genomes|Genetics. p. jkad217). While this approach has proven effective in improving optimization results, its application is constrained to optimizing only a limited number of parameters. As the number of parameters for optimization increases, the computational demands for performing sufficient simulations to obtain a broad coverage of the search space increases exponentially.

[0012] Recognizing these challenges, there is a need to create an optimization framework that requires fewer simulations without constraints in the number of optimized parameters and allows fully automation. A further issue for the optimization is that the evaluation of the objective function is computationally very expensive. Therefore the chosen optimization technique should allow several evaluations to be carried out in parallel.

[0013] Summary of the invention

[0014] The present invention meets this need by providing an enhanced optimization framework that builds on the concepts of kernel regression but additionally makes use of an evolutionary algorithm to allow for a more effective and general optimization. The methods of the invention consider a set of potential parameterizations of the breeding program, evaluate their performance based on stochastic simulations, and use these outputs to derive new parameterization to test in an iterative procedure. The evolutionary algorithm used achieves convergence to an optimum, typically the same optimum achieved in a conventional method (for example a prior art or known method or calculation), with a massively reduced number of simulations thereby obtaining much-reduced computing time, incorporation of class variables, and better scaling for a higher number of parameters considered in the optimization pipeline. This has been demonstrated herein exemplarily by comparison with the method described in Hassanpour et al. (supra), where 3 design parameters have been optimized.

[0015] The evolutionary algorithm framework has been designed to address the versatility of breeding programs with varying inputs capable of optimizing breeding programs with both continuous and discrete decision variables along the example of a dairy cattle breeding scheme (as also suggested in Hassanpour et al. (supra)), but is adaptable for any breeding scenario, regardless of the species, methodology, resources, or genetic traits involved — providing a versatile tool for a variety of breeding objectives and therefore is applicable to any plant or animal breeding program as long as it can be simulated / evaluated via stochastic simulations. A key aspect is that the described method allows optimizing a large number of parameters to enhance highly complex plant, such as wheat, breeding programs while maintaining equivalent budgets or penalizing higher costs by integrating the budget into the target function. The results reported herein demonstrate substantial improvements in genetic gain, with a 32.6% increase in genetic gain for the wheat line breeding program optimized solely for maximizing genetic gain. Even the diversity scenario for the wheat line breeding program, which balances genetic gain and genetic diversity for long-term sustain ability, achieved 4.5% more genetic gain than the baseline, while maintaining 9.1 % higher diversity. For the hybrid wheat breeding program, which optimizes more than 15 parameters, the optimized scenario increased genetic gain by 8.8% in the female component and 4.6% in the male component compared to the baseline. This makes the described approach particularly valuable for commercial breeders, as it enables efficient resource use, maximizes genetic improvement, balances between breeding goals, and enhances competitiveness, profitability, and sustainability.

[0016] In a first aspect, the present invention is thus directed to a method, in particular a computer-implemented method, for the optimization of a breeding program design, comprising the steps of:

[0017] (1) definition of the multivariate optimization problem and at least one termination criterion;

[0018] (2) initialization of the first set of parameterizations;

[0019] (3) evaluation of the first set of parameterizations based on the optimization problem defined in step (1);

[0020] (4) selection of a number of parameterizations from the first set of parameterizations based on the evaluation made in step (3) or from the second or further set of parameterizations based on the evaluation made in the preceding step (6), wherein said number of selected parameterizations is lower than the total number of parameterizations in the set of parameterizations;

[0021] (5) generation of a second or further set of parameterizations based on the selected parameterizations of step (4) comprising

[0022] (5a) creating new parameterizations from the selected set of parameterizations obtained in step (4) by combination of selected parameterizations; and / or (5b) creating new parameterizations from the selected set of parameterizations obtained in step (4) by modification of selected parameterizations;

[0023] (6) evaluation of the second or further set of parameterizations generated in step (5) based on the optimization problem defined in step (1) and derivation of optima, wherein unless the at least one termination criterion is met steps (4) to (6) are repeated; and

[0024] (7) assessment of the optima derived in step (6).

[0025] In various embodiments, step (1) comprises defining a search space for of the parameters to be optimized. The parameters to be optimized may comprise class variables and / or continuous variables. Step (1) may also comprise defining the breeding objectives (target function) and any potential constraints for the parameters to be optimized.

[0026] In various embodiments of the method, step (2) comprises creating initial parameterizations within the search space, for example using an initial population of possible parametrizations, and optionally also meeting potential constraints to the parameters to be optimized that have been previously defined. In various embodiments, steps (3) and (6) comprise simulating the respective breeding program with the selected parameterizations by use of a suitable simulation method, preferably by use of stochastic simulation.

[0027] In various embodiments, steps (4) to (6) are repeated multiple times.

[0028] Step (4) may, in various embodiments, comprise

[0029] (a) selecting parameterizations with the highest value of the objective function from the latest set of parameterizations;

[0030] (b) selecting parameterizations with the highest value of the objective function from a prior set of parameterizations; and / or

[0031] (c) selecting parameterizations from the area (of the search space) with the highest average value of the objective target function of the latest set of parameterizations.

[0032] In such embodiments, in step (4) sufficiently diverse parameterizations may be selected. Alternatively or additionally, step (4) may comprise at least alternative (c), preferably alternative (c) in combination with alternative (a) or alternative (b), more preferably may comprise all three alternatives (a) to (c). In various embodiments, for (c) a kernel regression method is used.

[0033] In various embodiments, step (5) can further comprise

[0034] (5c) selecting new parameterizations different from those in the latest set of parameterizations; and / or

[0035] (5d) selecting parameterizations that correspond to optima in prior iterations; and / or

[0036] (5e) selecting parameterizations that correspond to optima and promising settings (as determined, for example, by kernel regression) in the present iteration.

[0037] In various embodiments thereof, step (5) comprises at least alternatives (5a) and (5b), preferably (5a), (5b) and any one or more of (5c), (5d) and (5e).

[0038] In various embodiments, the derivation of optima in step (6) comprises using kernel density estimate to assess in which areas there are enough simulations to provide sufficient coverage and / or using kernel regression to estimate objective target function locally. In such embodiments, the derivation of optima in step (6) may comprise using kernel density estimate to assess in which areas there are enough simulations to provide sufficient coverage and using kernel regression to estimate objective target function locally in all other areas, wherein the parameterization with the highest value based on the kernel regressions is used as the optima.

[0039] In the methods described herein, the number of iterations of steps (4) to (6) may be determined by setting a threshold value for improvements of the objective target function by each iteration, wherein if said threshold value is not met, no further iteration is carried out.

[0040] In various embodiments, step (7) can comprise assessment of the suggested optima by simulating the suggested optima multiple times. The methods described herein may be used for design of a breeding program for plants, in particular crop plants, or animals, in particular livestock. Exemplary crop breeding programs in which the methods described herein may be applied include, without limitation, wheat, soy, cotton, and corn (maize) breeding programs. Such uses also form part of the invention.

[0041] In a further aspect, the invention is directed to a data processing system comprising means for carrying out at least steps (3) to (6), preferably steps (2) to (6), more preferably steps (1) to (7) of the method of the invention.

[0042] In another aspect, the invention also relates to a computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out at least steps (3) to (6), preferably steps (2) to (6), more preferably steps (1) to (7) of the method of the invention.

[0043] In still another aspect, the invention concerns a computer-readable data carrier having stored thereon the computer program of the invention.

[0044] The methods disclosed herein may be computer-implemented methods. Also encompassed by the present invention is a computer program product designed to perform the methods described herein.

[0045] Brief description of the drawings

[0046] Figure 1 provides a schematic summary of the overall pipeline of the developed evolutionary algorithm framework, representing key steps and their interconnections within the optimization process.

[0047] Figure 2 shows a dairy cattle breeding scheme, as used in the examples.

[0048] Figure 3 shows an example visualization of the Snakemake workflow. Shown are the rule names defined and their input-output relationships.

[0049] Figure 4 shows the frequency of expected outcome for Scenario 1 based on 100 replicates in (4a) for genetic standard deviations, (4b) for average kinship (Based on IBD).

[0050] Figure 5 shows the suggested optima for the individual parameters of the breeding program design for the number of test daughters (5a), test bulls (5b), and selected sires (5c), as well as a binary variable to control if additional budget is spent on phenotyping of test daughters (5d). The black horizontal line represents the estimated optima through comprehensive exploration, achieved by conducting over 100,000 simulations utilizing kernel regression (Hassanpour et al. (supra))

[0051] Figure 6 shows the performance of the suggested optima after each iteration assessed using kernel regression based on all simulations in Scenario 1 . The black horizontal line shows the estimated optima obtained from Hassanpour et al. (supra). Figure 7 shows the suggested optima across iterations for Scenario 2. Labels denote iteration numbers, while the line illustrates the iterative pathway. The dashed segment zooms in on the overlapping iterations within the optimal range from iteration 20 to 40. The black circle denotes the reference point for optimal parameter settings obtained through kernel regression, involving over 100,000 simulations.

[0052] Figure 8 shows the estimates of the optimum values in Scenario 2: (8a) Initial population, (8b) Iteration 5, (8c) Iteration 10, (8d) Iteration 15, (8e) Iteration 20, (8f) Iteration 25, (8g) Iteration 30, (8h) Iteration 35, and (8f) Iteration 40. The white points without border in the illustration represent simulations from more than five iterations back that are not considered as potential optima. White points with borders represent simulations that were excluded as optima due to a low kernel density estimation. Black points represent candidate optima from the kernel density estimation with the white point with a cross representing the finally chosen parameterization.

[0053] Figure 9 shows the share of binary variable to increase heritability of phenotyping being active in each iteration in (9a) for Scenario 3a with open symbols being the share of parents and filled symbols being the share of population, (9b) for Scenario 3b with open symbols being the share of parents (selected parameterizations in step (4)) and filled symbols being the share of population (parameterizations available for selection in step (4)).

[0054] Figure 10 shows schematically the initialization (step (2)) with a pool of parameter settings.

[0055] Figure 11 shows schematically the evaluation of the parameterizations based on the optimization problem / target function (step (3)).

[0056] Figure 12 shows schematically the selection based on high performance (step (4)).

[0057] Figures 13 and 14 show schematically the generation of new parameterizations (step (5)).

[0058] Figure 15 shows schematically the next iteration with new parameter settings (step (6)).

[0059] Figure 16 shows a schematic overview of a wheat line breeding program, highlighting key stages that utilize various cohorts: doubled haploid (DH), preliminary yield trial (PYT), advanced yield trial (AYT), and elite yield trial (EYT). The dotted line represents the selection of superior parents through recurrent genomic selection (GS), based on their genomic estimated breeding values (GEBV). The figure is taken from Bancic et al. (infra).

[0060] Figure 17 shows a schematic overview of the hybrid wheat breeding program. The pipeline illustrates the progression through key stages: doubled haploid (DH), headrows (HDRW), observational trials (OBS), testcross seed production (TC) and testcross yield trial (TC YT). Figure 18 shows the performance of the suggested optima for the wheat line breeding program after each iteration as assessed using kernel regression based on all simulations. The dark line indicates the EA gain breeding scheme and the light line represents the EA diversity breeding.

[0061] Figure 19 shows the suggested optima for the individual parameters of the wheat line breeding program design for number of A) crosses, B) DH lines produced per cross, C) entries per PYT, D) entries per AYT, E) entries per EYT, and F) new inbred parents to replace the oldest inbred parents. The horizontal line represents the reference scenario presented by Bancic et al. (infra). The dark lines indicate the EA gain breeding scheme and the light lines represent the EA diversity breeding.

[0062] Figure 20 shows the average increase in breeding values in genetic standard deviations across 100 replicates. The lowest band indicates the baseline breeding scheme, the highest band represents the EA gain breeding scheme, and middle band represents the EA diversity breeding scheme, all with the same budget.

[0063] Figure 21 shows genetic trends for wheat line breeding program with genetic standard deviation and. The middle line represents the baseline breeding scheme, the lower line represent the EA gain breeding scheme, and the upper line represent the EA diversity breeding scheme, all with the same budget.

[0064] Figure 22 shows the performance of the suggested optima for the hybrid wheat breeding program after each iteration as assessed using kernel regression based on all simulations.

[0065] Figure 23 shows suggested optima for the individual parameters of the hybrid wheat breeding program design with A) rep_TC2_YT f, B) n_Cross f, C) n_Crossm, D) n_DH f, E) n_DHm, F) n_OBS1f, G) n_OBS1m, H) n_OBS2f, I) n_OBS2m, J) n_OBS3f, K) n_OBS3m, L) n_OBS1f_Share, M) n_OBS2f_Share, N) n_OBS3f_Share, O) n_OBS2m_Share, P) n_OBS2m_Share, Q) n_OBS3m_Share. The lighter horizontal line represents the baseline breeding scheme, while the dark lines indicate the EA breeding scheme using the evolutionary algorithm.

[0066] Figure 24 shows the share of the binary variable to perform replication of the testcross yield trial on the female side rep_TC2_YTf in year 8 in each iteration for the optimized hybrid wheat breeding program.

[0067] Figure 25 shows the average increase in breeding values in genetic standard deviations across 100 replicates for A) n_Crossf and B) n_Crossm. Dark lines represent the EA breeding scheme using the EA, and lighter lines indicates the baseline breeding scheme with the same budget.

[0068] Figure 26 shows the level of heterozygosity based on the average of 50 replicates for A) n_Crossf and B) n_Crossm. Dark lines represent the EA breeding scheme using the EA, and lighter lines indicates the baseline breeding scheme with the same budget. Detailed description

[0069] Herein a comprehensive pipeline for optimizing breeding scheme designs using an evolutionary algorithm (EA) is provided. The EA is structured as an iterative process, wherein certain steps are reiterated until a termination criterion is met. For illustrative purposes, the individual steps of the algorithm are described herein using the same dairy cattle breeding program previously examined in Hassanpour et al. (supra) as well as wheat breeding programs. The invention is however not limited to such dairy cattle and wheat breeding programs, but is applicable to other animal or plant breeding programs.

[0070] “Breeding scheme design” or “breeding program design”, as used interchangeably herein, relate to the planned breeding of a group of animals or plants, usually involving at least several individuals and extending over several generations. Such programs or schemes typically are a collection of crossing, evaluation, and selection tasks and decisions which vary across breeding stages. As these are laborious and time-consuming, it is highly desirable to have tools to simulate and design such breeding strategies and apply the results thereof to the actual real-life breeding process.

[0071] As described herein above, the present invention is directed to a (computer-implemented) method for the optimization of a breeding program. Said optimization of a breeding program comprises at least the following steps (1) to (7):

[0072] (1) definition of the multivariate optimization problem and at least one termination criterion;

[0073] (2) initialization of the first set of parameterizations;

[0074] (3) evaluation of the first set of parameterizations based on the optimization problem defined in step (1);

[0075] (4) selection of a number of parameterizations from the first set of parameterizations based on the evaluation made in step (3) or from the second or further set of parameterizations based on the evaluation made in the preceding step (6), wherein said number of selected parameterizations is lower than the total number of parameterizations in the set of parameterizations;

[0076] (5) generation of a second or further set of parameterizations based on the selected parameterizations of step (4) comprising

[0077] (5a) creating new parameterizations from the selected set of parameterizations obtained in step (4) by combination of selected parameterizations; and / or

[0078] (5b) creating new parameterizations from the selected set of parameterizations obtained in step (4) by modification of selected parameterizations;

[0079] (6) evaluation of the second or further set of parameterizations generated in step (5) based on the optimization problem defined in step (1) and derivation of optima, wherein unless the at least one termination criterion is met steps (4) to (6) are repeated; and

[0080] (7) assessment of the optima derived in step (6).

[0081] Step (1) relates to the definition of the optimization problem. The optimization problem is related to a breeding program design problem or task. Accordingly, firstly, the breeding problem is formulated as an optimization problem. This can involve defining the breeding objectives or the objective target function, practical constraints, and (decision) variables that govern the breeding strategy.

[0082] The optimization problem may, for example, be defined as maximizing genetic gain or to balance genetic gain and diversity. The optimization problem may also be directed to certain traits, such as (grain) yield, recessive diseases (RD), and protein content.

[0083] Practical constraints are, for example and without limitation, budget constraints, number of crosses, spatial capacity, such as housing for cattle or field capacity for crops, number of parents replaced each year, etc.

[0084] The optimization problem is “multivariate” in that it comprises more than one variable to be optimized, typically 2 or more, 3 or more or 4 or more variables. Herein, it is demonstrated that the described method is particularly suitable and powerful to be used with multivariate problems that comprise 4 or more parameters, such as at least 5, at least 10, at least 15 or at least 20 variables or even more. These variables that are to be optimized are also referred to as “parameters” herein. Generally, in said step (1) two main types of decision variables may be considered. Firstly, variables with a limited number of possible realizations that can be viewed as categorical features of a breeding program. A non-limiting example of this is the question of whether genomic selection or marker-assisted selection is applied in a specific step of the breeding program. These variables are also referred to herein as class variables.

[0085] Secondly, variables that can take values from a continuous scale, or at least a large number of discrete realizations are considered. This may include the number of candidates selected, phenotyped, and genotyped, or weights in a selection index. For the sake of simplicity, these variables are referred to herein as continuous variables.

[0086] Step (1) can therefore include the definition of a search space for the parameters to be optimized, which may include the above-mentioned class variables and / or continuous variables. The search space may be limited by any potential constraints for the parameters to be optimized that have been defined or set by the user. The term “search space”, as used herein, relates to a defined range of values that all parameters to be optimized can take.

[0087] The terms “variables” and “parameters”, as used herein interchangeably, relate to features or values that can vary and influence the optimization problem, i.e. the breeding program design, as defined herein above. Exemplary variables include, without limitation, those used in the examples below, such as number of test progeny, number of test parents, number of crosses among parents to start a breeding cycle, number of DH lines produced per cross, number of entries per (preliminary yield / advanced yield / elite yield) trial, and number of new inbred patents selected each cycle based on genomic estimated breeding values (GEBVs) from the DH stage to replace the oldest inbred parents. The at least one termination criterion defines a condition that, once it is met, results in the method to be terminated. The termination criterion can, for example, be a preset number of iterations or may be defined via a threshold value for improvement of each iteration relative to the previous one that once it is no longer met, indicates that the method is to be terminated and / or the optimization can be considered complete. The termination criterion / criteria may be defined by user input and / or may be based on empirical data and / or an available database.

[0088] In various embodiments, step (1) of defining of the multivariate optimization problem and at least one termination criterion therefore comprises

[0089] (1) defining a search space for the parameters to be optimized;

[0090] (2) defining the breeding objectives (target function); and / or

[0091] (3) defining constraints for the parameters to be optimized.

[0092] In various embodiments, all three of the above are included in step (1).

[0093] Examples for such objectives, the constraints and the search space are given herein above and in the examples.

[0094] Step (2) is directed to the initialization of the first set of parameterizations. To initialize the evolutionary pipeline it is necessary to generate a starting population of potential breeding program designs to consider. In various embodiments of the method, step (2) comprises creating initial parameterizations within the search space, for example using an initial population of possible parameterizations, and optionally also meeting potential constraints to the parameters to be optimized that have been previously defined. This provides for a starting point for optimization.

[0095] “Parameterization”, as used herein, relates to defining and / or choosing parameters (variables), typically within a predefined search space that provides constraints to the parameters to be optimized. It also refers to a single set combination of the variables, i.e. an individual parameter setting where all variables / parameters of the optimization problem have been assigned a specific value. A “set of parameterizations”, as used herein, relates to the totality of all parameter settings used / simulated in this iteration cycle. This number of parameter settings, i.e. the number of settings within the set, is not particularly limited but may be selected by those skilled in the art based on the optimization problem.

[0096] In various embodiments, the initial parameterizations may be randomly generated within the search space or a pre-defined part thereof. For this, a generalized Bernoulli distribution for class variables and a uniform distribution for continuous variables may be used.

[0097] Depending on constraints, in some other embodiments, a different optionally non-random sampling procedure might be required, especially when the variables are interdependent, and logical dependencies for each decision may come from other decisions that share the same resources. In such embodiments, the user may define an initial set of parameterizations based on experience or other factors.

[0098] Irrespective on the methodology used, step (2) in various embodiments achieves the initialization of the first set of parameters by creating initial parameterizations within the search space and meeting the constraints to the parameters to be optimized. The creating may be done randomly or by a user, as described above.

[0099] In the next step (3) the suitability of each parameterization, meaning each breeding program design, is evaluated regarding is performance relative to the multivariate optimization problem, e.g. the breeding objective as given by the objective target function, typically by stochastic simulation.

[0100] For this evaluation, each respective breeding program (with the selected parameterizations ) can be simulated via stochastic simulation. Suitable simulation tools are known in the art and include, by way of example and without limitation, the R package MoBPS (Pook et al. (2020) Mobps - modular breeding program simulator. G3 Genes|Genomes|Genetics 10:1915-1918), AlphaSim (Faux et al. (2016) Alphasim: Software for breeding program simulation. The plant genome. 9), Adam (Liu et al. (2018) Adam-plant: A software for stochastic simulations of plant breeding from molecular to phenotypic level and from simple selection to complex speed breeding programs. Frontiers in plant science. 9:1926), and QMsim (Sargolzaei & Schenkel (2009) Qmsim: a large-scale genome simulator for livestock. Bioinformatics (Oxford, England) 25:680-681). It is however also possible to seamlessly integrate other simulation tools.

[0101] In various embodiments, to effectively address the computational challenge of the EA and particular step (3), the implementation of parallel processing and workflow management in an efficient manner may be crucial. For this purpose, automation tools that are known and available in the art may be used, for example the automation provided by the Snakemake workflow management system (Molder et al. (2021) Sustainable data analysis with snakemake. F1000Research. 10:33).

[0102] In step (4) parameterizations are selected. If step (4) is carried out for the first time, selection is made from the first set of parameterizations as (simulated and) evaluated in step (3). If step (4) is performed within an iterative cycle, i.e. is performed for the second or further time, selection is made from the parameterizations of most recent step (6) (with the exception that step (4b) allows selection from parameterizations from a previous iteration, i.e. not the most recent step (6) but an earlier step (6) or step (3)).

[0103] The general object of this step is to identify the most promising areas of the search space with respect to the optimization problem (i.e. the areas that provide for the best-performing parameterizations or those with the highest improvement over previous settings) to be further investigated in the following steps based on the results from the previous step. For this, typically simulated parameter settings that will be used as “parents” of the next iteration will be selected. The selection process can be made using different strategies. However, in various embodiments, the following three strategies are considered:

[0104] Highest value of the objective function (4a): The parent settings are selected based on the value of the objective function that was derived solely based on the simulation of the parameterization itself. This strategy ensures that the best performing parameterizations are prioritized in the process, thus enhancing the likelihood of generating successful, i.e. further optimized, parameterizations. The “highest value” as used in this context, refers to the best-performing parameterizations. “Best-performing”, as used in this context, refers to those parameterizations that give the best results with respect to the optimization problem. As typically more than one parameterization is selected, the term covers selecting a percentage or number of the best performing parameterizations. Said percentage may be, for example, 20 %, preferably 15% or 10%, more preferably 5% of the best results with respect to the optimization problem. The percentage may also be higher if desired, depending on the desired number or share of parameter settings to be selected in this step. If multiple parameterizations are selected, a diversity management strategy as detailed below may be used to avoid selecting parameterizations that are too similar to each other. The required dissimilarity is defined more concretely herein below and is similarly applicable to ensure that sufficiently diverse parameterizations from the set of parameterizations of the current iteration are selected. The term “sufficiently diverse” is used in this context to identify those parameterizations that have the predefined dissimilarity. Said dissimilarity may be preset using the Euclidian distance settings detailed below.

[0105] Highest value of the objective function from the previous iteration (4b): The parent settings are selected from previous iterations, insofar those have already been conducted. Typically, this means that best performing candidates from the last iteration are used. If, however, the expected performance of any suggested / predicted optima (see step (6)) is higher than all parameterizations of the current iteration, the most promising of these are added instead. This allows to reduce the risk of discarding potentially valuable solutions and high performing candidates due to stochasticity. Again, the “highest value” as used in this context, refers to the best-performing parameterizations. As typically more than one parameterization is selected, the term covers selecting a percentage or number of the best performing parameterizations from the given iteration.

[0106] Highest expected value of the objective function (4c): In this strategy, candidates based on the highest expected value of the objective function are selected. To determine this, kernel regression may be used, for example the method suggested in Hassanpour et al. (supra). This selection strategy has the advantage that instead of evaluating each candidate individually, kernel regression computes a weighted average of performance values for multiple candidates and is in contrast to other strategies not biased toward regions with more parameterizations tested overall. An adaptive bandwidth corresponding to the standard deviation in the individual parameter may be used herein. This means that parameterizations from the area with the highest average value of the objective function, as calculated using, for example, kernel regression, are used. The “highest average value” as used in this context, refers to the area that comprises, on average, the best-performing parameterizations. As typically more than one parameterization is selected, the term covers selecting a percentage or number of parameterizations within the area with the highest average value of the objective function. The term “area”, used in this context, relates to the part of the kernel regression computation with the highest average value. The selection then uses a (preset) number of parametrizations closest to the highest average value determined by kernel regression, optionally also using the predefined dissimilarity criterion. Typically there is more simulation of parameterizations with average values for the individual parameter, i.e. parameterizations that cluster in the center of the distribution. Therefore, the likelihood to get an outlier with respect to the value of the target function is higher in areas where less parameterizations are considered. Kernel regression helps with selecting settings from areas with low number of simulations and avoid early convergence to the center of the distribution. In various embodiments, step (4a) is carried out first and the kernel regression of step (4c) is applied to all values not selected in step (4a).

[0107] In addition, a diversity management strategy may be incorporated in the EA to avoid selecting highly similar parameter settings. In various embodiments, such dissimilarity criterion is defined as follows.

[0108] The iterative selection process assesses the Euclidean distance (d) between a candidate and the previously chosen parents. For each selection candidate, the distance to previously selected settings is calculated as follows: with: p being the number of parameters, and denoting the empirical standard deviation of the i-th parameter in the current iteration.

[0109] This condition ensures a candidate gets selected only if its minimum distance from the existing pool of settings is larger than a defined threshold. Generally, for (4b) the sole constraint may be that it has not been selected previously (>0). For step (4a) the threshold is typically lower than for step (4c), since the latter is less affected by stochastic variations and more similar settings tend to have more similar values smoothed objective function. As the variance in individual parameters decreases throughout iterations, it implicitly leads to the selection of more similar settings in later stages and a more focused evolutionary process.

[0110] In various embodiments, the step (4), at least strategy (4c) is used. In various preferred embodiments, at least two of the above or more preferably all three are used. This means that in various embodiments, alternative (4c) in combination with alternative (4a) or alternative (4b), more preferably all three alternatives (4a) to (4c) are used. These strategies may be weighted by the number of parameterizations selected accordingly. In various embodiments the number of parameterizations selected according to strategy (4a) is largest, followed by those selected according to (4c) and then (4b). In other embodiments, the number of parameterizations selected according to the used strategy is the same for each strategy. In one of the examples disclosed herein below, 30 parameterizations are selected, 20 of which are selected according to (4a), 7 according to (4c) and 3 according to (4b). Generally, the number of parameterizations selected may change over the iterations, for example in that in the first iteration a higher number is selected (as there may be more variation in the population) and said number may be (slowly) reduced over the following iterations. The ratio of the settings selected according to the different strategies may be chosen by those skilled in the art based on their routine knowledge and experience. In various embodiments, it is however preferred that at least strategies (4a) and (4c) are employed and the number of settings selected according to these strategies makes up at least 50% of the total number of selected settings, preferably at least 60, at least 70, at least 80 or at least 90% of the selected settings.

[0111] Accordingly, in various embodiments, in step (4) at least parameterizations from the area with the highest average value of the objective target function (4c) and / or with the highest value of the objective function (4a) are selected. The number of parameterizations selected by this approach may be in the range of at least 50% of the total number of parametrizations selected in this step. In various embodiments, the number may be higher, such as 60%, 70%, 80%, 90% or even 100%. In various embodiments, at least 20 % of all parameterizations selected are selected according to (4c), preferably at least 30, at least 40 or at least 50 %. In various embodiments, at least 20 % of all parameterizations selected are selected according to (4a), preferably at least 30, at least 40 or at least 50 %. In various embodiments, at least one of (4a) and (4c) is used, preferably both.

[0112] It is further generally preferred that the number of parameterizations selected relative to the total amount of simulated parameterizations is less than 50%, preferably less than 30 %, more preferably 20% or less, for example 19%, 18%, 17%, 16%, 15%, 14%, 13%, 12%, 11 %, 10%, 9 %, 8%, 7%, 6%, 5%, 4%, 3%, 2%, 1 % or less. Generally, it is understood that the number of parameter settings selected in step (4) is lower than the total number of simulations run and that there may be a minimum number of parameter settings that needs to be selected, for example at least 5 or at least 10 parameter settings.

[0113] In step (5) new parameterizations, i.e. the “offspring”, are created. For this, the previously selected parameterizations are used to generate a set of new parameterizations that are evaluated in the next iteration. This process involves applying various techniques to generate a new set of parameterizations. The number of created settings varies and can be set by those skilled in the art. The number may also vary based on which iteration the algorithm is in. For example, in the examples described herein below, in iterations 11 onwards a total of 300 settings are generated.

[0114] This step includes (5a) creating new parameterizations from the selected set of parameterizations obtained in step (4) by combination of selected parameterizations; and / or (5b) creating new parameterizations from the selected set of parameterizations obtained in step (4) by modification of selected parameterizations. These steps will be described in more detail herein below.

[0115] Initially, all previously chosen parameterizations may be considered again in the next iteration. This criterion facilitates assessing the same parameter settings, for example by using a new random seed for all parameterizations in the respective iteration, resulting in a more robust evaluation of their performance. The random seed or any alternative means for variation known to those skilled in the art is intended to result in a variation of the chosen parameterizations between different iterations even if the same parameterizations are used again, since the random seed / variation is different between the iterations. The random seed / variation may be determined by the user or may be based on empirical data or a preset mathematical function. Accordingly, it is generally possible to consider previously chosen parameterizations in this step a second or further time. However, at least part of the settings generated in at least one step (5) of the inventive method are new settings in that they differ from those of the current (and optionally all previous) iteration(s) and have been, for example, generated by step (5a) or (5b).

[0116] Step (5a) is directed to the combination of selected parameter settings. Specifically, in this step “new” parameterizations are created by combining multiple, typically two (optionally randomly chosen) parental parameterizations that have been selected in step (4). By using these parents, a new parameterization is created in each case. “New”, as used in this context means that the created parameterization is different from the parent parameterizations. In various embodiments, it may also be different from all previous parameterizations simulated.

[0117] For continuous variables this may be done by use of a weighted average, for example: wherein Xi and yi are the parent parameter settings and Zi is the newly generated parameter setting for parameter i. w may, in various embodiments, not be 0. The weighted average may be randomly selected.

[0118] For class variables, the newly generated parameter settings Zi are randomly sampled with an equal probability of belonging to either the class of xi or yi.

[0119] In some instances, it may further be necessary to round and scale continuous variables to obtain integer numbers while fulfilling constraints.

[0120] In later iterations of the algorithm, it can make sense to avoid combinations of parameterization with different values for the class variables as these settings might not be compatible with each other anymore, as non-class variables are adapted to work particularly well with the class variables.

[0121] In the present exemplified case, the probability of having two parameterizations with different class settings gets reduced as the algorithm progresses. By iteration 2, this probability is decreased by 20%, and it keeps decreasing by 10% per iteration until an 80% reduction is reached.

[0122] Step (5b) is directed to the introduction of small changes to selected parameters. This may be done alternatively or additionally to step (5a). In the example below, it is outlined for parameters that originate from the combination process in step (5a), but the claimed methods are not limited thereto. For continuous parameters, the size of the variation / in parameter / may be sampled from a uniform distribution with the range determined by the variance of the parameter and a mutation occurring in each respective parameter with a probability of pactiv - where oxi denotes the standard deviation and is derived based on the empirical variance of the parameterizations of the current iterations.

[0123] In the below described examples, in the first iteration, pactiv is set to 0.2, indicating a 20% chance of mutation for all parameters.

[0124] As detailed above, minor changes may be applied to the selected parent parameterization, i.e. not the combined new parameterizations, to generate new parameterizations according to step (5b). To avoid generating a setting already generated in the step of considering previously chosen parameterizations again, the mutation rates (pactiv ) may be increased (for example to 0.3) and the sampling procedure repeated in case no mutations are performed. If such step (5b), i.e. the modification of a parent parameterization is carried out, this is done with a probability of 100%, i.e. the above formula is modified to

[0125] Zmut = (zt + t-i, Z2 + 12, ...), for example with ~ U(-2aXi , 2oxi)

[0126] “Minor changes” or “small changes” as used herein relates to a range of possible modifications that is minor in that it does not drastically change the parameter but only slightly changes it. To achieve this, it may be determined by the (empirical) variance of the parameter / parameterizations of the current iteration. A measure for said empirical variance may be the standard deviation. The range of possible modifications may for example be set to be within double or triple the standard deviation range but may range over a broader spectrum such as 0.1 to 10 times the standard deviation, for example 0.2 times to 5 times or 0.5 to 5 times, typically 1 to 3 times, such as 2 times ((-2ox, , 2oxi).

[0127] In addition to strategies described as steps (5a) and (5b), which may be combined as described herein above, further options exist to generate new parameter settings that may also be comprised in step (5). These include, without limitation,

[0128] (5c) selecting new parameterizations different from those in the latest set of parameterizations; and / or (5d) selecting parameterizations that correspond to optima in prior iterations; and / or

[0129] (5e) selecting parameterizations that correspond to optima in the present iteration.

[0130] Step (5c) includes randomly selected new parameterizations within the search space. This is typically not done when calculating the site of mutations, i.e. when calculating based on the variance in settings. In embodiments that do not include such calculating the site of mutations, steps (5b) and (5c) may be jointly used. Steps (5d) and (5e) include selecting parameterizations that correspond to optima that have been derived in an earlier step (6), i.e. step (6) of an earlier iteration. Steps (5d) and (5e) may thus also be considered to use copies of the “parents”. It is however not preferred to include the same setting multiple times under the same random seed.

[0131] In various embodiments of step (5), said step comprises at least alternatives (5a) and (5b), optionally in combination in that combination and modification are performed on the same set of parameters, optionally in further combination with sets of parameters that are derived from either combination or modification. In various embodiments, step (5) comprises (5a) and (5b) and any one or more of (5c), (5d) and (5e).

[0132] Accordingly, in various embodiments, in step (5) parameterizations are generated by employing at least steps (5a) and (5b). The number of parameterizations selected by these two options (5a) and (5b) may be in the range of at least 50% of the total number of parametrizations generated in this step. In various embodiments, the number may be higher, such as 60%, 70%, or 80%. In various embodiments, at least 20 % of all parameterizations generated are generated according to (5a) and / or (5b), preferably at least 30, at least 40 or at least 50 %.

[0133] In step (6) the evaluation of the second or further set of parameterizations generated in step (5) based on the optimization problem defined in step (1) and a derivation of optima is carried out. Unless at least one termination criterion is met steps (4) to (6) are repeated. Steps (4) to (6) may be repeated multiple times, with the number of repeats either being limited by pre-defined constraints or the termination criteria, as defined below. In various embodiments, the number of repeats, also referred to herein as “iterations”, may range from 5 to 1000 or more, for example 10 to 500 or 10 to 300 or 10 to 100 or 15 to 80 or 20 to 60 or 20 to 50.

[0134] The optima may be derived by providing a graphical output of the optimized parameters over the course of the iterations and interpreting said graphical output either visually by the user or automatically by a computer program.

[0135] In various embodiments, to derive the optima, a kernel density estimation may be employed with the object to determine which areas of the search space include sufficient coverage. In some instances only those settings from the last iteration(s) with a value for the kernel density estimation above a given threshold level, such as of above 20% quantile of these parameterizations, may be considered to avoid using parameterizations in sparsely sampled areas. Generally, for the estimation of the kernel density estimation, only simulations from the last 2 to 10 iterations may be used, for example the last 5 iterations.

[0136] Exemplarily, this may be done using the following: with p being the number of parameters and using the empirical standard deviation in the last five iterations as the bandwidth hj and the use of a multivariate Gaussian kernel for K:

[0137] Subsequently, the expected performance of all remaining parameterizations may be estimated using a kernel regression (see Step 4c; (Hassanpour et al. (supra)) using all simulations. The parameterization with the highest value based on the kernel regression is then used as the optima.

[0138] In example 1 described herein below, 40 iterations were performed without any other termination criteria accessed, i.e. the termination criterion was set to be a preset number of iterations.

[0139] To define (different) termination criteria, the optima from all previous iterations based on a kernel regression based on all simulations can be assessed and if improvements on the target function and individual parameters for a given number of iterations, such as, for example, 5 or 10 or 15 iterations, are below a certain threshold terminate the pipeline. This threshold may be chosen depending on the specific optimization problem and can be highly dependent on the desired precision of results (Jain et al. (2001) On termination criteria of evolutionary algorithms. In Proceedings of the Genetic and Evolutionary Computation Conference: pp. 768-768; Ghoreishi et al. (2017) Termination criteria in evolutionary algorithms: A survey. In: Proceedings of 9th International Joint Conference on Computational Intelligence, Volume 1. pp. 40373-384. SCITEPRESS Digital Library).

[0140] Given the time-consuming nature of simulating real breeding programs, users may invest time in manually monitoring the algorithm’s performance, e.g. by visual inspection. If users observe no substantial changes or improvements over successive iterations in both parameter settings and the objective function’s value but also are considering the computational cost arising from the algorithm, the optimization process may be stopped earlier. “Substantial”, as used herein in relation to changes and improvements, means a change or improvement that is, depending on the optimization problem, considered to be high enough by the user to consider further iteration cycles for further optimization. The value for said improvement / change may be derived from a database or may be preset or predefined by the user or may be based on empirical data. A possible alternative approach could be to begin with a conservative iteration limit, i.e. set a number of iterations based on conservative estimates, and only increase it if the need arises after further evaluation. The objective is thus to repeat the assessment, selection, and generation of new parameter settings over successive iterations until a stable solution to the objective functions is achieved. Step (6) thus comprises, in various embodiments, evaluation of the second or further set of parameterizations generated in step (5) based on the optimization problem defined in step (1) and derivation of optima by simulating the respective breeding program with the selected parameterizations by use of stochastic simulation comprising using kernel density estimate to assess in which areas there are enough simulations to provide sufficient coverage and / or using kernel regression to estimate objective target function locally, wherein steps (4) to (6) are repeated until the at least one termination criterion is met. In some embodiments, both kernel density estimate to assess in which areas there are enough simulations to provide sufficient coverage and using kernel regression to estimate objective target function locally are used.

[0141] In step (7) the optima are assessed. This assessment involves a more detailed analysis / simulation of the optima determined by the preceding method steps. This analysis is necessary since the kernel regression carried out in step (6) will naturally be biased in an optimum. In this step (7), the suggested optimum may therefore be simulated a (high) number of times for further confirmation, for example 20 to 200 times. Step (7) thus serves as a verification tool to determine and confirm that the optima obtained in steps (1) to (6) provide for the desired results. Said verification may for example be done by defining a threshold by which the simulations carried out with the determined optimum values for the parameterization may vary. This accounts for potential variability in the outcomes and ensures robust estimates of key genetic descriptors, such as genetic gain and diversity.

[0142] Step (7) thus preferably includes simulating the suggested optima multiple times.

[0143] Only if this verification is successful, the results obtained may be applied and used in the actual breeding program. If this verification is unsuccessful, some additional iterations of the method may be carried out and / or the optima reassessed. In the (rare) event that the obtained values for the optimization are not accepted, the method may be repeated, typically with a different objective function and / or different constraints. The whole concept of the described method is schematically shown in the flowchart in Figure 1 .

[0144] Depending on the optimization problem and the iteration of the algorithm, adapting parameter settings can improve the convergence of the pipeline. The following provides some general guidelines on when and whether to deviate from the exemplified default and in which direction.

[0145] Regarding the initial population size (Step 2), the goal should be to obtain a good coverage of the initial search space. Therefore, with more parameters or larger search intervals, it can make sense to increase the size of the initial set of parameterizations from that exemplified herein below significantly, for example to a couple of thousand. In case simulations require a high computational load, it might be necessary to reduce the initial set of parameterizations. However, one should be aware that this will increase the risk of running into a local maximum.

[0146] To assess parameterizations more effectively, it may be advantageous to evaluate scenarios with multiple replicates within a single iteration (Step 3). This approach is particularly relevant for extremely small breeding programs and short time horizons where stochasticity can be a major factor in the evaluation. However, even for the comparably small example with three variables used herein below, such extensive replication was not necessary.

[0147] For selecting the best parameter settings (Step 4), it is recommended to select more parameterizations in the first couple of iterations, as there is still more diversity present and to avoid the loss of potentially promising settings. In the present example, 100 parameterizations have been selected in the first two iterations, 50 parameterizations in iterations 3-9, and 30 parameterizations afterward with similar splits between Steps (4a), (4b), and (4c) (see Table 1 below).

[0148] For generating new parameter settings (step 5), the same number of new parameterizations (in the present example: 300) in each iteration is created. However, as step 5.3 is mostly intended for fine- tuning already promising settings, this is not applied in the first few iterations, and more focus is given to steps 5.1 and 5.2. Regarding the more general disclosure of the inventive method, this means that in the first few iterations the focus is on using previous settings and using new settings that arise from combination of previous settings (and, optionally, additional modification).

[0149] “Promising settings”, as used herein, relates to those settings that provide desirable optimization results even if they are not (yet) optima. This may be determined by comparing them to previous settings wherein an improvement or an improvement beyond a preset threshold value indicates that these settings are promising.

[0150] As a refinement to Step 5, mutation / modification rates can be adjusted based on the observed changes in each parameter during previous iterations. For example, a binary parameter that consistently shows superior results with one of the parameterizations should undergo less frequent mutation.

[0151] In various embodiments, two or more parameters to be optimized may be linked in that if one parameter is modified there is a chance the a second ’’linked” parameter is changed as a result of the modification to the first parameter. This particularly applies to parameters that are modified (by mutation), as certain parameters are related to other parameters so that if one is modified the other(s) can or may need to be adapted, too. The method can thus be, in various embodiments, designed to account for such linked paramters in that a mutation / modification of one parameter has a chance (that can be set to a certain percentage) to cause an adaptation in at least one further parameter. Said chance can be set to any value but is preferably not set to 100% to avoid that the other parameters are always adapted. In various embodiments, said chance may be set to any value between 10 and 90%, for example between 20 and 80 %, such as 30, 40, 50, 60 or 70 %. In a specific embodiment, two to-be-optimized parameters may be linked such that if a mutation on one of the parameters to be modified occurs, there is a given chance, such as a 50% chance, that the other parameter is adapted to maintain the same total number of lines generated. To avoid inflation of the total number of ’’mutations” in these parameters, the ’’mutation rate” may additionally be reduced as desired, for example by 50%.

[0152] Further details are apparent from the more detailed description of step (5) in the examples below.

[0153] Disclosed herein is, inter alia, a novel EA framework developed to optimize breeding program design with both class and continuous design variables that are suitable for the joint optimization of multiple design parameters of breeding programs in a computationally efficient matter.

[0154] The methods described herein represent a much-enhanced version of previously reported pipelines, such as the kernel regression pipeline of Hassanpour et al. (supra). The iterative nature of the EA, where each iteration produces more data (simulation) for kernel regression allows for more reliable and efficient identification of suitable parameterization for a breeding scheme, with the kernel regression still used as a core element to cope with the challenges of an optimization problem with a target function that includes stochastics / noise in the evaluation of the target function as present when using stochastic simulations.

[0155] The developed methods provides a lot of flexibility to easily adapt parts of the algorithm to improve the efficiency of the algorithm but also the breeding program design itself. Nonetheless, models also provide robustness, as shown in exempla scenarios 1 , 2, and 3a with independent runs obtaining very similar final results. The term convergence should in this context be used with caution, as the finally obtained optima can gradually change with decreasing bandwidth. Therefore the suggested "optima" will most likely not be the exact optima but at least be very close to it. The stochasticity in the evaluation of the target function will naturally cause minor deviations between runs and although differences in solutions will exist, practically, suggested optima between the three scenarios that should have the same optima all performed very similarly. The suggested optima in Scenario 1 exactly matching the optima from Hassanpour et al. (supra) can therefore mostly be seen as a coincidence, but not as a sign of exact convergence.

[0156] Robustness is particularly highlighted by Scenario 2, demonstrating that the EA algorithm can find optimal solutions even outside of the initial search space. This is an advantage over the previously reported kernel regression method (Hassanpour et al. (supra)), which relies on predefined parameter bounds and cannot dynamically adapt its search space. As such, it falls short in terms of automation and efficiency. Naturally, the choice of suited parameters and design space in the EA will make convergence both quicker and more reliable.

[0157] Even Scenario 1 could have easily been improved in terms of required computations by not considering cases of more than 250 test bulls in the initialization. Nonetheless particularly in more complex optimization problems covering a wide range of the parameter space should usually be the priority. Particularly with a higher number of parameters, this should reduce the risk of running into local maxima. To avoid running into local optima, one could for example extend the generation of new parameterizations (step (5)) by randomly sampling parameterizations similar to the initialization. Most of the simulations generated this way will be highly explorative with a low likelihood of providing good solutions. Hence, greatly increasing computational load which at least in this relatively simple breeding scheme was not in support. Visual inspection of the change in target function and individual parameters is a common practice with evolutionary algorithms. In many cases, fine-tuning the parameters associated with termination criteria relies often on an iterative, trial-and-error approach (Jain et al. (supra)). By integrating Snakemake into the described EA framework, the flexibility in determining when to stop the algorithm was increased without the need to rerun initial iterations. It is understood that while Snakemake was used in this example, there are various alternatives known to those skilled in the art that can be used instead. It is however generally considered advantageous if the next simulation is run automatically whenever computing resources are available. In some instances, even if the value of the objective function remains relatively stable, there might be variation in individual parameter settings across iterations. This has practical implications in real world scenarios, and it enables breeders to potentially allocate resources differently or achieve the same outcome through alternative scenarios, which might be logistically more feasible.

[0158] In this regard, the definition of a suitable target function is of major importance, as from practical experience defining such an objective function is practically not easy or at minimum can very abstract. Here, visual manual human assessment can also help to check if the suggested optima from the EA is not only in the defined search space but also in the realm of solutions a breeder would grant reasonable I doable.

[0159] For practical breeding, it would also be conceivable to run the EA pipeline multiple times, potentially with different target functions to then perform an in-depth analysis to calculate key characteristics (genetic gain, inbreeding, etc.) from these optima and pick the preferred solution. By this, the abstract concept of a target function can be replaced with a practical choice between which combination of genetic gain I inbreeding is the most desirable. Different target functions could for example be generated by using different weightings of genetic diversity and short / long-term genetic gain.

[0160] Assessing the EAs’s performance in terms of speed and computational effort is broadly defined to include various sensible metrics, such as the number of iterations, CPU time, or any similar indicators, with the number of simulations required being the main driver of computational load in the described methods. The efficacy of the disclosed EA framework in reducing computational resources is highlighted through its performance across all scenarios. For example, in Scenario 1 only 2400 simulations were required to achieve similar results as compared to a previous kernel regression approach which relied on more than 100,000 simulations (Hassanpour et al. (supra)). This demonstrates a considerable reduction in computational effort, with the disclosed algorithm achieving a decrease of over 40-fold compared to the kernel regression method. This reduction holds even with a less-than-ideal initial search space or when adding a binary decision variable and therefore emphasizes the practical advantages of the disclosed EA-based method in optimizing large-scale breeding program designs.

[0161] Additionally, scaling for a higher number of parameters should be greatly improved as the traditional kernel regression scales exponentially in the number of assessed parameters.

[0162] In the context of economic considerations for optimization strategies, a crucial aspect involves evaluating the costs against potential benefits. For the described example, the total per-job cost is 0.55 cent per simulations. Running 40 iterations with 12.300 simulations would therefore result in a cost of 67.65€. In most commercial industrial practices, the focus of optimization lies in finding the best solution within a specific operating region or parameter space of interest that meets cost-effectiveness criteria and generates profits within a specified timeframe. In scenarios where a breeding program has limited prospects for improvement, allocating a substantial budget for optimization, may not be economically justified, particularly as the required computing time for larger breeding program simulation will be substantially higher.

[0163] In the process of optimizing breeding scheme design, Jannink et al. (Jannink et al. (2023) Insight into a two-part plant breeding scheme through bayesian optimization of budget allocations. Crop Science 2024: 1-12; DOI: 10.1002 / csc2.21124) showed that there are difficulties when using Bayesian optimization to allocate budgets effectively in breeding schemes. One of the limitations faced in this investigation and previous work (Hassanpour et al. (supra)) is the lack of support for class variables. The disclosed EA strategy helps address the computational intensity associated with continuous optimization problems involving class variables. These problems can be computationally intensive not only due to their combinatorial nature but also due to the increase in the number of possible outcomes. Particularly with a high number of class variables, the total number of combinations will rapidly increase knfor n class variables that all can take k values.

[0164] It is therefore recommended to use as low a number of class variables as possible, e.g., in the here considered Scenario 3b the additional cost of improved phenotyping was extremely low which from a human side makes it quite obvious to spend this additional money. However, as the resulting improvement for the target function is low, the overall upside is low, and high overall stochasticity in the evaluation is observed, the EA for all 40 iterations considered both binary settings. On the contrary, more substantial differences, such as an increase in the costs of phenotyping of 1000 Euro for a marginal improvement in precision, are more easily detectable by the algorithm to be unsuitable. Therefore, it is recommended that if such a design decision seems straightforward based on quantitative genetic theory or intuition, it might be advisable to simplify the optimization process by fixating such a class variable from the beginning or running separate optimization pipelines for both binary settings.

[0165] Furthermore, Jannink et al. (supra) reported high variability in the outcomes of different runs of the Bayesian optimization pipeline depending on small changes such as input genotypes and therefore lacking the ability to draw general conclusions from the obtained results. In the present methods, input genotypes and trait architecture were randomly sampled for each simulation. As very similar optima were obtained in the scenarios that should have the same optima (Scenarios 1 , 2, 3b) this should be a strong indicator of the generality and stability of the approach.

[0166] Described herein is thus an optimization framework using an EA that integrates a local search approach based on a kernel regression model which shows superior optimization efficiency to existing approaches and is applicable to both classes and continuous variables, thus enabling breeders to explore a wider range of scenarios compared to traditional methods.

[0167] The results across all problems indicate that the disclosed framework shows great promise by robustly estimated optima while significantly reducing computation time. It is further demonstrated that the EA algorithm consistently converges towards a common optimal solution, showcasing its robustness and ability to identify globally optimal or near-optimal configurations. The algorithm’s superior convergence speed, solution diversity, balance between exploitation and exploration, and robustness to stochasticity highlight its potential for larger breeding optimization tasks.

[0168] The adaptable nature of the disclosed framework makes it not only suitable for various future projects but also ensures flexibility in accommodating different breeding program designs. Users can easily modify, extend, or replace steps and adjust parameter choices as necessary. Thus, the disclosed framework supports optimization strategies that adjust to changing needs in breeding programs.

[0169] Examples

[0170] Example 1

[0171] Scenario 1 - Traditional dairy cattle scheme

[0172] All tests were executed on a server cluster with Intel Platinum 9242 (2X48 core 2.3 GHz) CPUs using Snakemake toolkit version 7.21 .0, which was configured to distribute single jobs via a SLURM scheduler to the backends of the cluster. Simulations were conducted on single nodes using a single core per simulation, taking approximately 15 minutes, and peak memory usage of 5 GB RAM per simulation. The computing time of all other steps combined increases approximately linearly in the number of iterations, but even in iteration 40 only took a negligible seven seconds.

[0173] Step (1):

[0174] A traditional dairy cattle breeding scheme, illustrated schematically in Figure 2 is considered. The three parameters that are considered for optimization are:

[0175] 1. xi: number of test daughter

[0176] 2. X2 : number of test bulls

[0177] 3. X3: number of selected sires For simplification purposes, explicitly not considered is the genotyping as commonly done in dairy cattle breeding in recent years. Further, only a single quantitative trait (milk yield, with heritability (h2) of 0.3) is considered.

[0178] As constraints, the breeding program at hand is limited by an annual budget of 10,000,000 EUR with housing costs of 3000 EUR per bull and 4000 EUR per cow. To avoid unnecessary computations, unreasonable settings like the use of too many bulls or selecting just a single bull per year are excluded by additional constraints: x1 + x2 + x3 > 0 100 < x2 < 700 3 < x3 < 30

[0179] 4000xi + 3000X2 - 10000000 < 0

[0180] Exactly as done in the previous study of Hassanpour et al. (supra), the objective function (m) is a linear combination of the expected genetic gain (g) and the expected inbreeding level ( f ) after 10 generations (5 years of burn-in + 10 years of future breeding), to prioritize / weigh between the genetic gain and diversity: m(x) = g(x) - 50 x f (x)

[0181] Subsequently, three additional examples are showcased to highlight the versatility of the developed EA framework, each representing a modification of the baseline (Scenario 1). These alternative scenarios are detailed below.

[0182] Step (2):

[0183] For the present toy example, an initial set of 600 parameterizations is generated with an additional scaling step on variables xi and X2 to ensure that the budget constraint is met. For instance, if xi = 1025, X2 = 300, and the total cost of the breeding program is 5,000,000 Euros. Given the budget is 10,000,000, with no advantages of underspending, both values are doubled. In case the budget constraint cannot be exactly met due to rounding, a breeding program with costs slightly below the maximum cost is used, by first rounding down and then marginally increasing parameters as much as possible.

[0184] Step (3):

[0185] The suitability of each parameterization selected in step (2) is evaluated regarding the breeding objective as given by the objective target function. For this evaluation, each respective breeding program was simulated via stochastic simulation using the R package MoBPS as described by Pook et al. (Pook et al. (2020). Mobps - modular breeding program simulator. G3 Genes|Genomes|Genetics 10:1915-1918).

[0186] Step (4): Simulated parameter settings are used as parents of the next iteration. The number of selected parents based on each step varies based on the iterations of the EA with values given below for iteration 11 onwards (see Table 1). In these iterations, 30 out of 300 settings are selected, as follows:

[0187] Step (4a):

[0188] 20 of the 30 of the parents are selected based on the value of the objective function that was derived solely based on the simulation of the parameterization itself. This guarantees that the best-performing parameterizations are prioritized in the reproduction process, subsequently enhancing the likelihood of generating successful parameterizations.

[0189] Step (4b):

[0190] 3 parameterizations from previous iterations are used. This typically means selecting the bestperforming candidates from the last iteration, however, if the expected performance of any suggested optima based on kernel regression (see Step (6)) is higher than all parameterizations of the current iteration, up to the three most promising of these are added instead. Hereby, the risk of discarding potentially valuable solutions and high-performing candidates due to stochasticity is reduced.

[0191] Step (4c): 7 candidates are selected based on the highest expected value of the objective function. For this, the kernel regression method suggested in Hassanpour et al. (supra) is employed. Such selection gives the advantage of instead of evaluating each candidate individually, kernel regression computes an expected value based on a weighted average of performance values for multiple candidates and in contrast to Step (4a) will not be biased toward regions with more parameterizations tested overall. In contrast to Hassanpour et al. (supra), here the use of an adaptive bandwidth corresponding to the standard deviation in the individual parameter is proposed.

[0192] Similar to real-world breeding programs, also incorporated is a diversity management strategy in the EA to avoid selecting highly similar parameter settings. Herein, to calculate the dissimilarity criterion, the iterative process assesses the Euclidean distance (d) between a candidate and the previously chosen parents. For each selection candidate, the distance to previously selected settings is calculated as follows: j - min

[0193] { k ( x> ) is a pre\ iously selected setting ) tm ' wherein p is the number of parameters, and

[0194] ' denotes the empirical standard deviation of the i-th parameter in the current iteration.

[0195] For step (4a), the threshold is 0.04 n, while given that step (4c) is less affected by stochastic variations and more similar settings tend to have more similar values smoothed objective function, the threshold is raised to 0.08 n. In contrast, for step (4b), the sole constraint is that it has not been previously selected (> 0). As the variance in individual parameters decreases throughout iterations, it implicitly leads to the selection of more similar settings in later stages and a more focused evolutionary process.

[0196] Table 1

[0197] Note: Step (4a) presents individual simulations with the highest objective function value from the current iteration, Step (4c) presents individual simulations with the highest expected value of the objective function determined by kernel regression, and Step (4b) presents individual simulations with the highest objective function value from the previous iteration and previous optima. Step 5.1 represents the sum of Steps (4a), (4b), and (4c). Step 5.2 represents the number of new settings (simulations) generated through a combination of selected parameterizations. Step 5.3 represents the number of new parameterizations created through minor modifications of selected settings.

[0198] Step (5): Generation of new parameterizations

[0199] Subsequently, the previously selected parameterizations are used to generate a set of new parameterizations that are evaluated in the next iteration. This process involves applying various techniques to create a new set of parameterizations, drawing inspiration from the process of meiosis, as described in the following. The number of generated settings will again vary based on which iteration the algorithm is in (see Table 1). The values given below are for iterations 11 onwards with a total of 300 settings generated.

[0200] Step (5.1): Selected parameter settings

[0201] Initially, all 30 previously chosen parameterizations are considered again in the next iteration. This criterion facilitates assessing the same parameter settings using a new random seed in each iteration, resulting in a more robust evaluation of their performance.

[0202] Step (5.2): Combination of selected parameter settings

[0203] Furthermore, 180 parameterizations are generated by combining two randomly chosen parental parameterizations in Step (4), denoted as = (xi, X2, ...) and Y = (yi, y , ...). By using these parents, a new parameterization Z = (zi, Z2, ...) is created in each case.

[0204] For continuous variables this is done by the use of a weighted average: z, = wxi + (1 - w)y, with w ~ U(0, 1)

[0205] For class variables, z, is randomly sampled with an equal probability of belonging to either the class of Xi or yi. In this example, it is necessary to furthermore round and scale continuous variables to obtain integer numbers while fulfilling the budget constraint (see Step (1)). Following the combination process, small changes are introduced to the combined parameterizations, inspired by the process of allelic mutation in meiosis.

[0206] For continuous parameters, the size of the variation f, in parameter / is sampled from a uniform distribution with the range determined by the variance of the parameter and a mutation occurring in each respective parameter with a probability of pactiv - and mactiv.i ~ B(0, Pactiv ) where oxi denotes the standard deviation and is derived based on the empirical variance of the parameterizations of the current iterations.

[0207] In the used algorithm, in the first iteration, pactiv is set to 0.2, indicating a 20% chance of mutation for all parameters. In later iterations of the algorithm, it can make sense to avoid combinations of parameterization with different values for the class variables as these settings might not be compatible with each other anymore as non-class variables are adapted to work particularly well with the class variables. In the present case, the probability of having two parameterizations with different class settings gets reduced as the algorithm progresses. By iteration 2, this probability is decreased by 20%, and it keeps decreasing by 10% per iteration until an 80% reduction is reached. Similarly, the mutation rate pactiv in class variables is reduced by 10% in iteration 2 and reduced by 5% per iteration until a reduction of 40% is reached.

[0208] Step (5.3): Minor modifications of selected parameter setting

[0209] Lastly, the selected parent parameterizations are considered and minor changes / mutations applied to them (see Step 5.2). To avoid generating a setting already generated in Step 5.1 , mutation rates pactiv are increased to 0.3 and the sampling procedure is repeated in case no mutations are performed.

[0210] To further improve convergence speed there are various minor improvements to consider regarding mutation rates. For class variables, in each iteration, the share of each class is calculated for both the pool of parameterizations and the selected parents. If for the last five iterations, the sum of the share of all but a single class is lower than the highest mutation rate (Step 5.3: pactiv = 0.3) and the share of these classes in the selected parameterizations in Step 5 is lower than 0.75 * pactivj = 0.225, mutation rates for this parameter are halved.

[0211] If conditions are still fulfilled for the reduced mutation rate, the mutation rate is further reduced to 0.01 .

[0212] For continuous variables, additional prioritization and direction are given to some parameters instead of symmetrically sampling mutations around prioritizing one direction: with psign on default being 0.5. To prioritize, the kernel regression is used to approximate the expected change in case of an increase or decrease in the parameter:

[0213] In case either of the two is positive, pSignis increased / decreased by 0.03 to make it more likely to sample in that respective direction.

[0214] In case both nr, and m'n are negative, the overall mutation probability m, is reduced by 20%. This procedure is repeated for the optima from the last five iterations.

[0215] In case at least five continuous parameters are considered, the local derivates for all parameters in the current iteration are compared and the mutation rate mi in the 20% of the parameters with the most favorable mutation rate are increased by 50%. Contrarily, the least favorable mutation rates are reduced by 50%. In case a mutation rate exceeds 0.3 / 0.4 for Step 5.2 / 5.3 it is reduced to this respective threshold.

[0216] To avoid non-impactful mutation in case a parameter is fixed, is also advised to specify a minimum value for the mutation range that is then used instead of 2 sigma,. For all scenarios in this study, a minimum range of 40, 20, and 3 have been used for xi, X2, and X3, respectively.

[0217] Step (6): Convergence / Optima / Termination criteria

[0218] To derive the optima, first a kernel density estimation was employed to determine which areas of the search space include sufficient coverage. Herein, it is proposed to only consider those settings from the last five iterations with a value for the kernel density estimation above the 20% quantile of these parameterizations, to avoid using parameterizations in sparsely sampled areas. For the estimation of the kernel density estimation, only simulations from the last five iterations are used. with p being the number of parameters and using the empirical standard deviation in the last five iterations as the bandwidth hj and the use of a multivariate Gaussian kernel for K:

[0219] K(yi, yP) = Ki(yi) ... Kp(xp) with

[0220] Subsequently, the expected performance of all remaining parameterizations was estimated using a kernel regression (see Step 4c; (Hassanpour et al. (supra)) using all simulations. The parameterization with the highest value based on the kernel regression was then used as the optima.

[0221] In this example, 40 iterations were performed without any termination criteria accessed, i.e the termination criterion was the number of iterations. To define a general termination criteria, it is proposed to assess the optima from all previous iterations based on a kernel regression based on all simulations and if improvements for ten iterations are below a certain threshold terminate the pipeline. This threshold needs to be chosen depending on the specific optimization problem and is highly dependent on the desired precision of results (Jain et al. (supra); Ghoreishi et al. (supra)). Given the time-consuming nature of simulating real breeding programs, it is suggested that users invest sufficient time in manually monitoring the algorithm’s performance, e.g. by visual inspection. If users observe no substantial changes or improvements over successive iterations in both parameter settings and the objective function’s value but also the computational cost arising from the algorithm, they may consider stopping the optimization process earlier. A possible alternative approach could be, to begin with a more conservative iteration limit and only increase it if the need arises after further evaluation.

[0222] Step (7): Final assessment of the optima

[0223] After the optimum is identified, the finally obtained breeding scheme is analyzed in-depth, as the kernel regression will naturally be biased in an optimum. For this, the suggested optimum is simulated a high number of times, e.g. in the present case 100 replicates.

[0224] Alternative Scenarios

[0225] Scenario 2 - Reduced initial search space

[0226] In a previous study (Hassanpour et al. (supra)), optimal parameter settings (2368, 175, 19) for the three parameters also considered herein in Scenario 1 were identified. To assess the effectiveness of the EA algorithm disclosed herein, Scenario 2 tackles the same resource allocation problem as Scenario 1 .

[0227] However, an initial search space that does not contain the optima was considered to determine the ability of the EA to still identify the optima. For this the following two additional constraints were considered in the initialization (Step 2):

[0228] 300 < X2 s 500

[0229] 15 < X3< 25

[0230] Scenario 3 For showcasing the optimization of a class variable, Scenario 1 was extended by introducing a binary variable, X4. This variable represents a breeding strategy that leads to improved phenotyping that leads to a reduced residual variance and hence a higher heritability (h2= 0.32). For this study, it was left open as to what this new breeding strategy is, but one could envision more uniform housing conditions, the use of electronic devices to measure physiological status, or large-scale collection of additional data like mid-infrared spectroscopy. Step 3 is accordingly adapted by sampling initial values forx4 from a Bernoulli distribution X4 ~ B(0.5). Here two versions of the scenario (3a / 3b) are considered with varying additional costs of 1 ,000€ / 10€, respectively. Resulting in the new constraints:

[0231] 3a : xi(4000 + 1000x4) + 3000x2- 10000000 < 0

[0232] 3b : xi(4000 + 10x4) + 3000x2- 10000000 < 0

[0233] Snakemake

[0234] The EA algorithm disclosed herein is iterative, involving multiple interdependent steps, where the computational demands for executing some steps are notably high, particularly the resource-intensive simulation of breeding programs in Step 3. To effectively address this computational challenge, the implementation of parallel processing in an efficient manner is crucial. For this purpose, the present optimization pipeline makes use of the automation provided by the Snakemake workflow management system (Molder et al. (2021) Sustainable data analysis with snakemake. F1000Research. 10:33).

[0235] An illustrative representation of the Snakemake workflow for the described evolutionary optimization model is provided in Figure 3.

[0236] The Snakemake process is built around four rules, corresponding to the steps of the algorithm, which are initialization (step 2), evaluation (step 3), evolutionary algorithm (steps 4, 5, and 6), and the final in- depth analysis of the obtain optima (step 7).

[0237] To clarify, the individual simulations within Step 3 are completely independent of each other and can easily be run in parallel using the built-in capabilities of Snakemake to seamlessly interact with various job scheduling systems (e.g., SLURM https: / / slurm.schedmd.com / documentation.html) on distributed hardware stacks, thus ensuring portability to a wide range of hardware setups. For configuring this setup, reference is made to the Snakemake plugin catalog at https: / / snakemake.github.io / snakemake-plugin- catalog / index.html. This catalog provides a starting point for configuring Snakemake to work with various cluster schedulers, ensuring optimal distribution and execution of multiple tasks across the computing cluster.

[0238] Results

[0239] Application of the described evolutionary pipeline to the optimization problem formulated in Scenario 1 suggests a final optimum of 2368 test daughters and 175 test bulls, of which 19 test bulls are selected (x1 = 2368, x2 = 175, x3 = 19) with an expected outcome for the target function of 107.041 , with a genetic gain of 9.07 genetic standard deviations (Figure 4a) and increase of the inbreeding level of 0.0426 (Figure 4b) after 10 generations based on the averages of 100 replicates of this scenario.

[0240] All three individual parameters very quickly reached values close to the finally suggested optima (Figure 5a-5c). When evaluating the suggested optima per iteration based on all simulations conducted (to avoid effects of reduced bandwidth over iteration) even after seven iterations a value of 107.038 for the target function based on kernel regression was obtained (Figure 6), despite kernel regression by design being downward biased in the optimum.

[0241] In Scenario 2, basically the same optima of xi = 2374, X2 = 167, X3 = 19 after 40 iterations is obtained. Even though, the suggested optima in the first couple of iterations are performing slightly worse, the value for the target function is on par with the optima suggested in Scenario 1 after seven iterations (Figure 6). For the individual parameters, more change can be observed with the initial iterations suggesting using more test bulls, due to limitations in the initial search space. Similarly, a higher number of sires are selected (X3) to use a similar selection intensity. In the first ten iterations of Scenario 2, the primary emphasis is on quickly bringing X2 into its optimal range due to the higher overall impact of X2 on the optimum. Overall, more change in the individual parameters is observed with X3 in iteration 10 being as low as 15 (Figure 5c). A detailed overview of the changes in X2 and X3 is given in Figure 7.

[0242] Although parameterizations with smaller X2 values do exist in early iterations (Figure 8b), these are not considered in deriving the optima (Step (6)) due to the low kernel density. In later iterations the area of solutions considered as the optima more and more shifts towards the area with the expected optima (dark grey area, Figure 8c:8i).

[0243] In Scenario 3a basically the same optima as in Scenarios 1 & 2 is obtained with the binary not being active (x = (2365, 179, 20, 0)). In terms of convergence speed, more variation in individual parameters is observed in early iterations, particularly with iterations 4 and 5 suggesting optima that are later identified as poor solutions (Figures 5a:5c & 6) caused by the stochasticity in the evaluation of the target function (Step 3).

[0244] In contrast, the binary is active in Scenario 3b which allows for a higher overall value of the target function of 107.200 that represents a statistically significant improvement based on a t-test (Student 1908) (p < 0.00663). Due to the higher overall housing costs, the number of cows and bulls is slightly decreased with a finally suggested optimum of x = (2361 , 175, 20, 1). Regarding the binary parameter, mutation rates in Scenario 3a are reduced to half from iteration 17 onwards, while for Scenario 3b mutation rates remain high for the entire 40 iterations.

[0245] In both scenarios, the share of the favorable binary in the selected parameterizations (Step 3) is higher than in the overall population but both cases are considered over the entirety of the simulations with 18% / 25% of the parameterizations of the respective alternative binary setting (Figures 9a and 9b). Example 2: Optimization of a wheat line breeding program

[0246] In this study and example 3 below, two distinct wheat breeding programs using two different stochastic simulators were evaluated. State-of-the-art breeding schemes were selected as baseline scenarios and serve as effective benchmarks for assessing the performance of the inventive methods.

[0247] A wheat line breeding scheme has been used before by Bancic et al. (Crop Science 2024. doi:10.1002 / csc2.21312) to illustrate the advantages of integrating genomic selection (GS) early in the breeding process serves here as a baseline for existing breeding practices. AlphaSimR, a stochastic simulator (supra) was employed to evaluate and compare various breeding strategies. The simulations focused on two distinct target functions: one solely optimizing genetic gain in yield, and another balancing both genetic gain and genetic diversity. A summary of this baseline breeding program, its sizes, and costs are given in Figure 16 and Table 2 below.

[0248] In the baseline breeding scheme, each breeding cycle begins with initial crosses in the first year from 50 potential parents. In the second year, the focus is on generating F1 / doubled haploid (DH). In the second year, rapid recurrent GS is applied to DH lines, and parents are recycled for the next cycle. From the third to the sixth year, the program advances selected lines to yield trials, involving growing the plants in different environments to evaluate their performance and identify the best-performing lines. At the end of each breeding cycle, the best-performing line will be released as a new variety. Each cycle will collect the information of all the advanced lines from PYT (preliminary yield trial), AYT (advanced yield trial), and EYT (elite yield trial) to train the GS model with two years of data. In the given baseline breeding scheme, 20% of inbred parents from a total number of 50 parents in each breeding cycle were replaced using newly selected parents, chosen based on their highest genomic estimated breeding values (GEBVs) at the DH stage. The total cost associated with this breeding scheme amounts to 503,500 USD in each breeding cycle.

[0249] Table 2. Cost of the baseline wheat line breeding scheme based on Bancic et al.. For abbreviations, please refer to the description of Figure 16. For the optimization of two breeding goals in the wheat line breeding scheme, the evolutionary algorithm used in Example 1 has been adapted to the above wheat line breeding scheme. In the following subsections, the individual steps of this adaptation for the specific use case are describe in more detail.

[0250] Step 1 : Definition of the optimization problem

[0251] In this example, focus was on optimizing six variables with

[0252] 1 . n_Cross, the number of crosses among parents to start a breeding cycle

[0253] 2. n_DH, the number of DH lines produced per cross

[0254] 3. n_PYT, the number of entries per preliminary yield trial

[0255] 4. n_AYT, the number of entries per advanced yield trial

[0256] 5. n_EYT, the number of entries per elite yield trial

[0257] 6. n_ParentsReplace, the number of new inbred parents selected each cycle based on genomic estimated breeding values (GEBVs) from the DH stage to replace the oldest inbred parents.

[0258] Given an annual budget limit of 503,500$ and considering the expenses associated with each stage (Table 2), the practical constraints are presented as follows: n_Cross, n_DH, n_PYT, n_AYT, n_EYT, n_ParentsReplace > 0 n_Cross > n_DH n_PYT > n_AYT > n_EYT n_Cross < 500 n_PYT < 1000 n_ParentsReplace < 50 n_Cross x 60 + n_Cross x n_DH x 45 + n_PYT x n_AYT x 675 + n_EYT x 1000 < 503,500

[0259] For logistical and practical reasons, the number of n_Cross is set to a maximum of 500 lines (e.g., costs associated with making more crosses, managing more field plots, handling additional labels and packets, collecting more data, and conducting more parental checks (Witcombe J, Virk D. Euphytica. 2001 ;122:451-62. doi:10.1023 / A:1017524122821). At the n_PYT stage, the maximum field capacity is limited to 1000 lines. The number of parents replacing each year (n_ParentsReplace) cannot exceed the total number of the 50 oldest parents.

[0260] Two distinct breeding objectives were evaluated for the line wheat breeding program. The first strategy's objective function aimed to maximize the expected genetic mean of DH lines by year 20, prioritizing only genetic gain. This scenario is referred to as "EA gain breeding scheme" herein. rn = max E | / (x, ^)] gain (xGX)

[0261] The second strategy balanced genetic gain and diversity using a composite target function. This scenario is referred to as "EA diversity breeding scheme" herein. m = max E [0.8 * g(x, ) + 0.2 * (20 * f (x, ())] diversity (xGX)

[0262] In the above, m is the set of six optimal decision variables for the optimization problem, X denotes the feasible set of solutions for the decision variables, E represents the expectation, which is the result of the stochastic simulation over a realization of the random variable Xi, g is the expected mean genetic gain of DH, and f is the genetic standard deviation of DH lines after 20 years of breeding, as an indicator for genetic variance in the population. This function assigned 80% weight to mean genetic gain and 20% to maintaining genetic standard deviation. Factor of 20 is applied for appropriate scaling, ensuring that the genetic mean component is comparable in magnitude to the genetic standard deviation component.

[0263] Step 2: Initialize the first set of parameter settings

[0264] The decision variables have been subjected to predefined ranges during the initial sampling of relevant parameter settings. All variables were drawn from a uniform distribution (Table 3). Variables may exceed or fall below their initial boundaries during the optimization process if they are not constrained by hard limits in Step 1 . After sampling the first set of breeding program designs, all continuous variables should be scaled to ensure they exactly match with the budget constraint, as outlined above. In this step, an initial set of 1000 parameter settings (breeding program designs) is generated.

[0265] Table 3. Initial range of decision variables for optimizing the wheat line breeding program. For abbreviations, please refer to the description of Figure 16.

[0266] Step 3: Evaluate new settings

[0267] To assess the effectiveness of different breeding program designs, stochastic simulation was used to model each program and evaluate its performance to the breeding objective. The simulation script was taken from Bancic et al. (supra) (available at https: / / github.com / HighlanderLab / jbancic_alphasimr_plants) and adapted to the EA framework as suggested above.

[0268] Step 4: Select parameter setting

[0269] The same methodology was employed for selecting optimal parameter settings during the optimization process as outlined in Example 1 . Reference is also made to Table 4 below. Table 4. Parameter setup for the number of selected and generated settings for wheat line and hybrid wheat breeding scheme step designation is according to Figure 1).

[0270] In the following the designation of steps is in accorndance with Figure 1. Step 4.1 presents selection of individual simulations with the highest objective function value from the current iteration, Step 4.2 presents selection of individual simulations with the highest expected value of the objective function determined by kernel regression, and Step 4.3 presents selection of individual simulations with the highest objective function value from the previous iteration and previous optima. Step 5.1 represents the sum of Steps 4.1 , 4.2, and 4.3. Step 5.2 represents the number of generated new settings through a combination of selected parameter settings. Step 5.3 represents the number of new parameter settings created through minor modifications of selected settings. Replication indicates the number of replications per simulation for each iteration.

[0271] Step 5: Generate new parameter settings

[0272] The approach for generating new parameter settings during the optimization process mirrored the methodology established in Example 1 (See Table 4 above). The total number of lines at the end of year two in the breeding scheme results from n_Cross x n_DH instead of a single parameter. To enable ’’mutations” to be weighted between n_Cross and n_DH used without modifying the total number of lines generated, these two parameters were modeled as linked parameters. This means that if a mutation on one of the parameters occurs, there is a 50% chance that the other parameter is adapted to maintain the same total number of lines generated. To avoid inflation of the total number of ’’mutations” in these parameters, ’’mutation rates” were additionally reduced by 50%. This provides an extension to the original concepts in Exampe 1 (Step 5.3), where only ’’mutations” to individual parameters were considered.

[0273] Step 6: Stabilization / Optima / Termination criteria

[0274] Due to the fast computation time of the simulations, the EA framework for the wheat line breeding scheme was configured to continue for 150 iterations. No specific termination criteria were established.

[0275] Step 7: Final assessment of the optima

[0276] The final assessment is executed, once the EA reaches its final iteration. The resulting EA breeding scheme, after identifying the ideal values for each optimized parameter, is thoroughly analyzed, and the results are compared with the baseline breeding scheme. For this comparison, the expected genetic mean, genetic standard deviation, and selection accuracy at the DH stage per breeding cycle with constrained costs over 20 years is evaluated based on data from 100 independent runs for each scenario. Results

[0277] Applying the evolutionary framework to the optimization problem to maximize two target function in the wheat line breeding program (Figure 1), the following optimal parameters were identified after 150 iterations (Figure 18; Table 5).

[0278] Table 5: Suggested optima for the individual parameters of the wheat line breeding program design in iteration 150. The best value represents the baseline scenario presented by Bancic et al. (supra). For abbreviations, please refer to Figure 3.

[0279] The optimization algorithm exhibited a clear convergence pattern, with the objective function stabilizing after iteration 50 in the EA gain breeding scheme and after iteration 90 in the EA diversity breeding scheme (Figure 19). Despite fluctuations in individual parameters (Figure 18), further iterations did not result in substantial improvements, indicating that an optimal solution was achieved.

[0280] The number of initial crosses increased by 301 % to 401 in the EA gain scheme and by 175% to 275 in the EA diversity scheme, compared to the baseline of 100 crosses (Figure 19:A). The number of DH lines per cross decreased in both optimized schemes, which is of primary interest, with the EA gain scheme showing a more dramatic reduction of 84.27% to 14 lines, while the EA diversity scheme decreased by 64.04% to 32 lines, compared to the baseline of 89 lines (Figure 19:B). This resulted in a total of 5,614 DH lines for the EA gain scheme (a 36.92% reduction) and 8,800 for the EA diversity scheme (a slight 1.12% decrease) from the baseline's 8,900 lines.

[0281] The allocation of lines across yield trials also differed significantly. In the PYT, the EA gain scheme decreased the number of lines by 18.6% to 407, while the EA diversity scheme increased it by 35.4% to 677, compared to the baseline of 500 (Figure 19:C). For AYT, the EA gain scheme showed a substantial increase of 360% to 230 lines, while the EA diversity scheme decreased by 48% to 26 lines, from the baseline of 50 (Figure 19:D). In EYT, the EA gain scheme increased by 40% to 14 lines, while the EA diversity scheme decreased by 50% to 5 lines, from the baseline of 10 lines (Figure 19:E). The proportion of parents replaced per cycle also varied greatly, with the EA gain scheme increasing by 360% to 46 (92%), while the EA diversity scheme maintained the baseline level of 10 (20%) (Figure 19:F). After identifying the suitable parameters, the results are presented as an average across 100 simulation replicates (thick lines) also with 95% confidence interval (light bands) for each breeding scheme demonstrate significant improvements in the optimized breeding scheme compared to the baseline (Figure 20). The wheat line breeding program optimized solely for maximizing genetic gain demonstrated a 32.6% increase in genetic gain compared to the baseline scheme, demonstrating the effectiveness of resource reallocation with the same budget. The positive difference in the genetic mean is due to the EA gain scheme emphasizing greater genetic diversity through increased crosses, coupled with a substantial increase in parent replacement, which indicates intensified selection pressure. The EA diversity focused scenario achieved a 4.5% higher genetic gain than the baseline, even with a more balanced approach, demonstrating that an optimized strategy, including diversity, outperforms the nonoptimized baseline breeding scheme.

[0282] The use of genomic selection with more emphasis on maximizing genetic gain has led to a more rapid decline in genetic variance across all scenarios (Figure 21 : A), with the selection accuracy for both EA breeding schemes and the baseline breeding scheme being relatively similar (Figure 21 : B). Regarding genetic variance (Figure 21 : A), by year 20, the baseline breeding scheme maintained genetic variance with a genetic standard deviation of 0.353, representing a 49.6% decrease from the initial value in year 1. The EA gain breeding scheme showed a more pronounced reduction, with genetic variance expresssed in standard deviation decreasing by 73.9% to 0.183. In contrast, the EA diversity breeding scheme slightly increased genetic varinace expressed in standard deviation to 0.385, representing a 9.1 % increase, compared to the baseline at year 20.

[0283] Optimization and computing time for simulation

[0284] Simulations for the wheat line breeding program, conducted using AlphaSimR simulator Version 1.5.3, required approximately 1 minute and 0.5 GB of peak memory usage per simulation on a single core. One replication per simulation was performed during the optimization process. The EA framework was conducted using the Snakemake workflow management system (version 7.21.0), which distributed individual simulations through a SLURM scheduler to the cluster backends.

[0285] Example 3: Optimization of a hybrid wheat breeding program

[0286] Structure of baseline breeding program

[0287] The hybrid wheat breeding program illustrated in Figure 17 is inspired by the benefits of dividing breeding objectives into two separate parts, as outlined by Gaynor et al. (Crop Science. 2017;57:2372- 86. doi:10.2135 / cropsci2016.09.0742). The population improvement component aims to quickly enhance the population's average genetic quality through repeated GS. The product development component is focused on identifying and developing new hybrid varieties, ensuring that the improved genetic material from the population improvement component is used to create market-ready varieties but is ignored in the exemplary simulations (i.e. no parameter was used to make this part variable). The following section describes the baseline breeding scheme, with the numbers provided in each year corresponding to this scenario outlined in Table 6. Table 6

[0288] Y = year; St. = stage; fem. = female; pop. = population; Env. = environments; Sei. = selection; Gen. = genotype cost / unit (USD); Phen. = phenotype cost / unit (USD).

[0289] Cost and key features of the hybrid wheat baseline breeding scheme with doubled haploid (DH) lines, headrow trials (HDRW), observational trials (OBS), testcross seed production (TC), testcross yield trials (TC YT), replication of testcross yield trials (rep TC2 YT). The considered traits are grain yield (GY), recessive disease (RD), and protein content (PC). The heritability h2values represent estimated approximate values based on breeding value estimation, considering the fluctuating residual variance observed across various stages of the breeding program. * denotes that 15 hybrids are produced in total, drawing from the identified 30 best female and male DHs.

[0290] Simulation of the hybrid wheat breeding program Base population and burn-in phase The breeding program's structure is based on winter wheat Triticum aestivum L., as described by Gaynor et al. (supra). The genome was simulated for each of the 21 chromosomes, each with a genetic length of 1 .43 Morgans. The resulting genome sequences had 1 ,000 bi-allelic SNPs per chromosome (21 ,000 total). Simulated genome sequences were used to generate 200 founder lines, partitioned into two separate heterotic pools. In order to obtain our base population, the hybrid wheat breeding program begins with the filling process once the breeding program stages are defined and the founder population is simulated (also termed burn-in, see Gaynor et al.).

[0291] Each breeding stage is filled with a unique cohort that traces back to its parent lines and replicates the overlapping breeding cycles. For this stage, within each pool, two homozygous founder testers were also simulated in the initial stage of the program.

[0292] This filling process mimics the simultaneous operation of different breeding cycles in real-world programs, which means filling each stage of the breeding pipeline with a distinct cohort, all originating from the same parental population but at different stages of development. In the hybrid wheat breeding program with five stages: 1) Cross, 2) DH lines, 3) observational trials 1 (OBS1), 4) observational trials 2 (OBS2), 5) and observational trials 3 (OBS3), five distinct cohorts are required. The filling process involves progressing the first cohort through all stages, ultimately saving it in the final stage (OBS3), which becomes the oldest cohort, with subsequent cohorts following at earlier stages of the pipeline.

[0293] To effectively simulate a more realistic representation of genetic gain across generations, reflecting the continuous efforts to improve traits in breeding populations, the filling process ensures that later generations are progressively better than their predecessors. These cohorts are simulated in a way that avoids a scenario where generations 1 to 5 have essentially the same genetic values. A detailed description of the simulation procedure in MoBPS can be found in Chapter 9.15 of the MoBPS User Manual (MoBPS supra).

[0294] Three traits were simulated, assuming each trait was controlled by 300 underlying purely additive and 30 dominant QTLs, with effect sizes drawn from a Gaussian distribution. Trait 1 can be seen as a grain yield (GY), trait 2 can be seen as a protein content (PC), and trait 3 represents a recessive disease (RD). Traits were assigned weights according to their importance for the desired breeding aims. GY has an index weight of three units per genetic standard deviation (gSD), making it the most important trait, while PC and RD both have an index weight of 1 unit per gSD. QTL effects were standardized to have an underlying genetic mean of 100 and a genetic variance of 10.

[0295] Following Pook et al. (bioRxiv. 2025:2025.01.10.632416. doi:10.1101 / 2025.01.10.632416) traits were simulated with separate realizations in male and female parental lines (per se performance) and in hybrids (cross performance). The correlations between per se and cross performances were set at 0.2 for PC and 0.1 for RD. GY was considered only as a hybrid trait (cross performance). The fully considered correlation matrix between traits is shown in Table 7 and a detailed explanation of correlated traits can be found in Chapter 15.1 of the MoBPS (supra). Furthermore, the phenotyping of individual traits in 10 different environments was considered. Genotype-by-environment interactions are simulated by assessing different variants for the different traits depending on the environment. Therefore, for each trait by environment, a separate trait is simulated, resulting in 50 trait combinations when using 5 traits in 10 environments. The required correlation matrix (50x50) is generated by combining the correlation matrix between traitstrait(5x5) and environmentsCTlv(10xl0) using the Kronecker Product:

[0296] 2 = Strait ® 2 env

[0297] As no real data for environments was available in this study, the correlation matrix between environments was previously sampled by sampling pairwise correlations for a predefined range (0.6 - 0.8). As pairwise sampling does not necessarily lead to a positive semi-definite matrix, the resulting matrix was subsequently projected in the space of positive semi-definite matrices by setting all negative eigenvalues to 0 (matrix. posdef in MoBPS).

[0298] Table 7. Expanded correlation matrix of traits and the relationships between parental lines (per se) and hybrid (cross) levels for grain yield (GY), protein content (PC), and recessive disease (RD), including their interactions.

[0299] PC (Per se) RD (Per se) GY (Cross) PC (Cross) RD (Cross)

[0300] PC (Per RD (Per se) GY (Cross) PC RD (Cross) se) (Cross)

[0301] For the simulation of the marker-assisted selection (MAS) process in MoBPS, a subset of QTLs associated with the trait of interest are randomly selected, and their effects on the trait are estimated. Predicted phenotypic values are then calculated using these genetic markers and their estimated effects, by selecting a fraction of the total QTLs and using linear regression methods to achieve a realistic prediction accuracy in this stage. To ensure that the heritability (h2) values accurately reflect the varying residual variance across different stages of the breeding program, simulated phenotypes were generated directly with the appropriate values presented in Table 3. A detailed explanation can be found in Chapter 5.4.1 of the MoBPS (supra).

[0302] At the start, 100 crosses each are performed in the female and male parental pool every year (Table 6). Parental combinations in each pool are chosen randomly from all possible combinations for the 100 parental lines in the crossing block. 100 DHs are produced per cross.

[0303] In year 4, DH lines are planted in headrows, where breeders visually assess plant development, including diseases, and discard undesirable lines based on field observations. Lines are also genotyped to apply MAS, ensuring early fixation of key traits at an early stage. Assessments of MAS along with visual agronomic evaluation are made in one location. 500 lines are selected after disease evaluation and MAS and enter replicated, multi-location observation trials (OBS1 ; Table 6).

[0304] OBS1 presents the evaluation of DH lines across two environments for the RD trait, while genomic prediction is simultaneously applied to all OBS1 individuals for GY and PC. The training population for GS included all available hybrid data of the past three seasons. The training population was continuously updated in subsequent years by incorporating new yield trial evaluations as they became available. At the same time, the same 500 DH lines / pool are engaged in testcross productions with 2 testers from the opposite pool to furnish seeds for the subsequent testcross yield trials (TC1 , Table 6). The 300 best DH lines are selected based on the per se performance of RD trait and their predicted general combining ability (GCA) for GY and PC traits to advance lines to the first testcross yield trials (TC1 YT, Table 6).

[0305] 300 best DH lines / pool are next tested again in observational trials (OBS2) for RD and their derived hybrids are tested in testcross hybrid trials for GY, PC, and RD (TC1 YT, Table 6). Concurrently, the same 300 DH lines / pool are engaged in testcross productions with 2 testers from the opposite pool to furnish seeds for the subsequent testcross yield trials (TC2, Table 6). The top 30 parents / pool are selected based on their observed and estimated GCA for all traits within a multi-environment trial (derived from TC1 YT). Additionally, selection considers per se performance for RD, derived from OBS2. The top 30 DH lines / pool are tested further in observational trials (OBS3) for RD and PC, and their derived hybrids are tested in testcross hybrid trials for GY, PC, and RD (TC2 YT, Table 6). Finally, hybrid are produced making selected combinations with the 30 best females and 30 best males (hybrid prod, Table 6). The top 15 parents / pool are selected on GCA for all traits (derived from TC2 YT) and additionally on the per se performance for RD and PC (derived from OBS3).

[0306] Retesting female hybrids

[0307] In the female breeding pool, a subset of testcross hybrids undergoes additional evaluation in subsequent testcross yield trials (repTC2 YT, Table 6) to enhance the accuracy of EBVs.

[0308] Recycling of parents

[0309] In accordance with Pook et al (supra), parents are recycled from multiple stages. 20% of parents consist of new materials that could, e.g., come from the characterization of new materials in the pre-breeding program. 60% of parents are derived from the OBS3 stage (twice evaluated for TC hybrid yield), while the remaining 20% is evenly split, with 10% coming from OBS2 (once evaluated for TC hybrid yield) and another 10% from OBS1 (only predicted GCA). In this selection process, top-performing candidates with highest EBVs were allowed to be used approximately five times more often for crosses than the lowest- ranked lines.

[0310] Objectives and parameters for optimization

[0311] Step 1 : Definition of the optimization problem To assess the overall effectiveness of the hybrid wheat breeding program, 17 parameter designs were selected for optimization as follows:

[0312] 1. rep_TC2_YT f whether to perform replication of the testcross yield trial on the female side in year 8

[0313] 2. n_Crossf being the crosses on the female side to start a breeding cycle

[0314] 3. n_Crossmbeing the crosses on the male side to start a breeding cycle

[0315] 4. n_DHf being the total number of DH lines produced on the female side

[0316] 5. n_DHmbeing the total number of DH lines produced on the male side

[0317] 6. n_OBS1f as the number of entries in observational trial OBS1 produced on the female side

[0318] 7. n_OBS1mas the number of entries in observational trial OBS1 produced on the male side

[0319] 8. n_OBS2f as the number of entries in observational trial OBS2 produced on the female side

[0320] 9. n_OBS2mas the number of entries in observational trial OBS2 produced on the male side

[0321] 10. n_OBS3f as the number of entries in observational trial OBS3 produced on the female side

[0322] 11 . n_OBS3mas the number of entries in observational trial OBS3 produced on the male side

[0323] 12. n_OBS1f_Share as the proportion of recycling parents selected from the OBS1 female cohort

[0324] 13. n_OBS2f_Share as the proportion of recycling parents selected from the OBS2 female cohort

[0325] 14. n_OBS3f_Share as the proportion of recycling parents selected from the OBS3 female cohort

[0326] 15. n_OBS1m_Share as the proportion of recycling parents selected from the OBS1 male cohort

[0327] 16. n_OBS2m_Share as the proportion of recycling parents selected from the OBS2 male cohort

[0328] 17. n_OBS3m_Share as the proportion of recycling parents selected from the OBS3 male cohort

[0329] Considering the 2,363,550$ budget limit for the entire breeding scheme and the expenses linked to each stage (Table 6), the following practical constraints are outlined:

[0330] All continuous variables > 0 n_OBS3f> 15 & n_OBS3m> 15 n_OBS1f_Share + n_OBS2f_Share + n_OBS3f_Share = 0.8 x n_Crossf n_OBS1m_Share + n_OBS2m_Share + n_OBS3m_Share = 0.8 x n_Crossm(n_Crossf + n_Crossm) x 30 + (n_DHf + n_DHm) x 80 + (n_OBS1f + n_OBS1m)x620 + (n_OBS2f+ n_OBS2m) x 240 + n_OBS3mx490 + n_OBS3fx (490 + 50 x rep_TC2_YTf) < 2,363,550

[0331] In hybrid wheat breeding, it is crucial to test a sufficient number of lines to ensure reliable yield performance and stability estimates in multi-environment trials before final release. For this practical reason, the number of n_OBS3f and n_OBS3mis set to a minimum of 15 lines and cannot go below this constraint. It is understood that this is an arbitrary choice and may be different in other simulations. The shares of recycled parents selected from the OBS1 , OBS2, and OBS3 cohorts were constrained and standardized to sum up to 80% of the total number of crosses on the female and male sides each, ensuring that 20% of new genetic material is introduced each year. In the simulation of the hybrid wheat breeding program, the focus was on the core genetic aspects and omitted certain practical steps to streamline the model. For this, seed production stages were not included in the simulation, as these do not directly impact genetic calculations. Similarly, while actual hybrid production was not simulated, its associated costs were incorporated into the overall cost function.

[0332] For this study, an objective function was established to maximize the average total genetic gain achieved over the 20 years of breeding. In many commercial breeding programs, immediate returns are more valuable than future returns. This idea is captured by applying an interest rate (r) to the hybrid wheat breeding program, similar to the concept in economics. In this study, r was set to 0.05 (5%) in the objective function, meaning that genetic progress made sooner is weighted more heavily than progress made later:

[0333] Where AG is the average total genetic gain over 20 years between the male and female part of the breeding scheme. The genetic gain for male and female populations is calculated using the following formulas: where A Gfemaieand AGmalerepresent the genetic gain achieved in the female and male populations in year t, based on the true underlying genomic values of the individuals, and r is the interest rate, set to 0.05.

[0334] Step 2: Initialize the first set of parameter settings

[0335] The initial bound of decision variables during the sampling process is as detailed below with all continuous variables drawn from a uniform distribution and the binary variable drawn from a Bernoulli distribution (Table 8). Similar to the wheat line breeding program, scaling was applied to ensure that all selected breeding program designs remained within a defined budget. Regarding the initial population size for the hybrid wheat breeding program, which involves a larger number of parameters, the goal is to achieve good coverage of the initial search space. Therefore, an initial set of 2000 parameter settings was created.

[0336] Table 8. Initial range of decision variables for optimizing the hybrid wheat breeding program. For abbreviations, refer to Table 6.

[0337] Step 3: Evaluate new settings

[0338] Breeding program designs were assessed using the MoBPS stochastic simulator with a simulation script available at https: / / github.com / AHassanpour88 / Evolutionary_Snakemake in the same git of the EA folder with exemplary scripts.

[0339] Step 4: Select parameter setting

[0340] The number of selected parameter settings can be found in Table 4.

[0341] Step 5: Generate new parameter settings

[0342] The process of generating new parameter settings for the hybrid wheat breeding program was consistent with the approach used for the wheat line breeding program, as detailed in Table 4. To facilitate ’’mutations”, a linked parameter was implemented to ensure that the combined shares of recycled parents from the 0BS1 , 0BS2, and 0BS3 cohorts always amounted to 80% of the total number of crosses. To prevent an increase in the total number of ’’mutations” in linked parameters, the ’’mutation” rates for the linked parameters were further reduced by 50%.

[0343] Step 6: Stabilization / Optima / Termination criteria

[0344] Given the high computational cost per simulation and the large number of parameters to optimize, the EA framework for the hybrid wheat breeding scheme was set to run for a maximum of 100 iterations, with no specific termination criteria other than a visual assessment.

[0345] Step 7: Final assessment of the optima

[0346] Following identifying the optimal solution through visual analysis of the proposed optima, the resulting EA breeding scheme was subjected to a detailed analysis, and its outcomes were compared to those of the baseline breeding scheme. To evaluate the effectiveness of the inventive framework, the mean genetic gain for the n_Crossf and n_Crossmwas determined per breeding cycle by averaging their true breeding values. Similarly, to show genetic diversity in each population, the share of heterozygous markers for the n_Crossf and n_Crossmwere evaluated over breeding cycles between scenarios. Optimization and computing time for simulation

[0347] Simulations for the hybrid wheat breeding program, using MoBPS Version 1.11.64, required about 40 minutes and 8 GB of peak RAM usage per simulation when run on 2 cores. One replication per simulation was performed during the optimization process. The EA framework was conducted using the Snakemake workflow management system (version 7.21.0), which distributed individual simulations through a SLURM scheduler to the cluster backends.

[0348] Results

[0349] The impact of optimizing resource allocation in a complex hybrid breeding program was investigated, considering both the female and male sides while prioritizing short-term genetic gain. Applying the described EA framework led to the final optimal parameters in Table 9. Although the termination criterion was set to 100 iterations, the algorithm was stopped at iteration 75 due visual inspection showing no change in the objective function value. Despite the program’s complexity with 17 parameters, the objective function stabilized after iteration 30 (Figure 22) and remained unchanged from iteration 70 onward, with all parameters reaching stable values (Figure 23).

[0350] Table 9: Suggested optima for the individual parameters of the hybrid wheat breeding program design in iteration 75. For abbreviations, refer to Table 6.

[0351] The optimization process led to substantial changes in resource allocation between the baseline and the optimized hybrid wheat breeding scheme. Notably, the n_Crossf decreased significantly by 76%, from 100 to 24 (Figure 23:B), whereas the n_Crossmexperienced a moderate reduction of 13%, from 100 to 87 (Figure 23:C). For the total number of DH lines, while n_DHf increased slightly by 11.05%, from 10,000 to 11 ,052 (Figure 23:D), the n_DHmdeclined by 5.3%, from 10,000 to 9,469 (Figure 23:E).

[0352] Adjustments were also observed in the allocation of observation plots across different stages. The number of first-stage female observation plots (n_OBS1f) saw a minor rise of 6.6%, from 500 to 533 (Figure 23:F), while the male counterpart (n_OBS1m) decreased by 10.4%, from 500 to 448 (Figure 23:G). In the second stage, the n_OBS2f remained relatively stable, decreasing slightly by 2%, from 300 to 294 (Figure 23:H), while the n_OBS2mdeclined by 32.3%, from 300 to 203 (Figure 23:l). In the third stage, n_OBS3f were reduced by half, from 30 to 15 (Figure 23:J), whereas n_OBS3mshowed a marginal increase of 3.3%, from 30 to 31 (Figure 23:K).

[0353] The share of first-stage female observation plots (n_OBS1f_Share) decreased from 60% to 42% (Figure 23:L), while the second-stage share (n_OBS2f_Share) more than doubled, rising from 10% to 21 % (Figure 23:M). Similarly, the third-stage share (n_OBS3f_Share) increased from 10% to 17% (Figure 23:N). On the male side, n_OBS1m_Share declined slightly from 60% to 55% (Figure 23:0), while n_OBS2m_Share remained unchanged at 10% (Figure 23:P), and the n_OBS3m_Share saw a modest rise from 10% to 15% (Figure 23:Q).

[0354] Additionally, in the baseline breeding scheme of the hybrid wheat breeding program, an additional replication of rep_TC2_YTf raises the costs associated with that step in the eighth year of the breeding program on the female side. However, it is important to note that the overall budget remained fixed in that changes were compensated for by scaling cohort sizes and only this specific step incurred more expenses. The best strategy for the selected optima was identified as choosing not to perform rep_TC2_YTf (Figure 23:A), with no changes observed from iteration 14 onward.

[0355] Figure 24 illustrates that during the optimization process, certain iterations favored a specific simulation with the binary setting activated, temporarily increasing its representation in the population. However, this preference was not sustained, as subsequent iterations of EA did not select this setting, causing its representation to decline. Nevertheless, across all simulations, 1 % of the generated parameter settings still reflect this alternative binary setting (Figure 24).

[0356] Expected genetic gain (gSD units) per cycle on the female side is primarily observed in generation 1 , due to a reduction in the generation interval (Figures 25:A). In contrast, gains on the male side demonstrate more consistency over time (Figures 25:B). The EA breeding scheme achieved 0.45 more gSD for n_Crossf (8.8%) and 0.24 more gSD for n_Crossm(4.5%) in year 20 compared to the baseline scenario (Figure 25). The comparison between the EA and baseline breeding schemes showed significant differences for both n_Crossf and n_Crossm.

[0357] The genetic gain per step for each scenario, based on n_Crossf and n_Crossmwithin one breeding cycle, is illustrated in Table 10, along with the corresponding selection intensities. The EA scheme consistently achieved higher genetic gain (genetic standard deviations; gSD) than the baseline breeding scheme across all steps for the female side. When considering consecutive steps for n_Crossf, the cumulative genetic gain in the baseline scheme from OBS2 to OBS3.1 and OBS3.1 to OBS3.2 was 1 .26 gSD, which remained slightly higher than the single-step gain of 1.22 gSD achieved by the optimized scenario. However, the small additional gain from OBS3.1 to OBS3.2 in the baseline scheme was compensated by a more efficient allocation of selection pressure in earlier stages, allowing the EA scheme to maintain a higher overall genetic gain. On the male side, the most substantial difference was observed in the OBS1 to OBS2 step, where the optimized scenario in n_Crossmapplied a much higher selection intensity (32.7% vs. 60% in the baseline) and achieved 0.14 gSD more genetic gain.

[0358] Table 10. Comparison of selection intensity and the amount of genetic gain for each step of the breeding programs in both the baseline and EA breeding schemes. EA<f> referes to optimization results on the female side, where EA<m) presents the results on the male side.

[0359] Over 20 years, both the baseline and EA breeding schemes led to a slight reduction in diversity on both the female and male sides. By Year 20, the baseline scheme retained slightly more diversity indicated by the share of heterozygosity (0.298 for females and 0.297 for males) compared to the EA breeding scheme (0.283 for females and 0.284 for males). This suggests that while the optimized breeding scheme achieved higher genetic gain, it also resulted in a marginally faster decline in genetic diversity. However, the small differences indicate that the integration of external genetic material (20%) helped sustain diversity to some extent in both approaches.

[0360] In examples 2 and 3, two optimized scenarios obtained from the EA framework with different decision variables were compared against a baseline breeding scheme, maintaining equivalent budgets for wheat line and hybrid breeding programs. The results allowed for a direct comparison of breeding strategies in terms of expected genetic gain at the same investment level. The results demonstrate the effectiveness of the proposed EA framework and provide valuable insights for its application in complex breeding designs.

[0361] The examples demonstrate the successful application of a stochastic simulation and evolutionary optimization framework to three distinct, complex breeding programs. A key strength of the claimed method is that it is the first to apply an advanced optimization framework to a complex breeding program, optimizing over 15 interdependent design parameters, including both continuous and class variables. The framework’s adaptability is evident through its application across various line and hybrid programs, using different breeding simulators. It is shown that, with a modest investment in implementing the present framework, breeding programs can achieve significant improvements in genetic gain. This approach offers considerable value for modern breeding programs aiming to accelerate genetic gain and maximize return on investment by refining existing strategies that require resource reallocation. The findings suggest that optimizing natural trade-offs, such as allocating resources for phenotyping versus genotyping and evaluating a larger number of candidates versus conducting more thorough assessments of fewer candidates, can substantially improve efficiency without additional financial investment. However, if desired, financial adjustments can be made to either reduce or increase the cost of the program, providing flexibility for further optimization and improvement of breeding outcomes.

[0362] Finally, the use of the EA algorithm enables more strategic planning, allowing breeders to explore scenarios, predict outcomes, and adjust their strategies as needed.

Claims

Claims1. Computer-implemented method for the optimization of a breeding program design, comprising the steps of:(1) definition of the multivariate optimization problem and at least one termination criterion;(2) initialization of the first set of parameterizations;(3) evaluation of the first set of parameterizations based on the optimization problem defined in step (1);(4) selection of a number of parameterizations from the first set of parameterizations based on the evaluation made in step (3) or from the second or further set of parameterizations based on the evaluation made in the preceding step (6), wherein said number of selected parameterizations is lower than the total number of parameterizations in the set of parameterizations;(5) generation of a second or further set of parameterizations based on the selected parameterizations of step (4) comprising(5a) creating new parameterizations from the selected set of parameterizations obtained in step (4) by combination of selected parameterizations; and / or(5b) creating new parameterizations from the selected set of parameterizations obtained in step (4) by modification of selected parameterizations;(6) evaluation of the second or further set of parameterizations generated in step (5) based on the optimization problem defined in step (1) and derivation of optima, wherein unless the at least one termination criterion is met steps (4) to (6) are repeated; and(7) assessment of the optima derived in step (6).

2. The method of claim 1 , wherein step (1) comprises defining a search space for the parameters to be optimized, wherein the parameters to be optimized preferably comprise class variables and / or continuous variables, wherein step (1) preferably further comprises defining the breeding objectives (target function) and constraints for the parameters to be optimized.

3. The method of claim 1 or 2, wherein(i) step (2) comprises creating initial parameterizations within the search space and meeting potential constraints to the parameters to be optimized; and / or(ii) steps (3) and (6) comprise simulating the respective breeding program with the selected parameterizations by use of a suitable simulation method, preferably by use of stochastic simulation; and / or(iii) steps (4) to (6) are repeated multiple times.

4. The method of any one of claims 1 to 3, wherein step (4) comprises(a) selecting parameterizations with the highest value of the objective function from the latest set of parameterizations;(b) selecting parameterizations with the highest value of the objective function from a prior set of parameterizations; and / or(c) selecting parameterizations from the area with the highest average value of the objective target function of the latest set of parameterizations.

5. The method of claim 4, wherein(1) in step (4) sufficiently diverse parameterizations are selected; and / or(2) step (4) comprises at least alternative (c), preferably alternative (c) in combination with alternative (a) or alternative (b), more preferably comprises all three alternatives (a) to (c).

6. The method of claim 4 or 5, wherein for (c) selecting parameterizations from the area with the highest average value of the objective target function of the latest set of parameterizations a kernel regression method is used.

7. The method of any one of claims 1 to 6, wherein step (5) further comprises(5c) selecting new parameterizations different from those in the latest set of parameterizations; and / or(5d) selecting parameterizations that correspond to optima in prior iterations; and / or(5e) selecting parameterizations that correspond to optima or promising settings, preferably determined by kernel regression, in the present iteration.

8. The method of claim 7, wherein step (5) comprises at least alternatives (5a) and (5b), preferably (5a), (5b) and any one or more of (5c), (5d) and (5e).

9. The method of any one of claims 1-8, wherein the derivation of optima in step (6) comprises using kernel density estimate to assess in which areas there are enough simulations to provide sufficient coverage and / or using kernel regression to estimate objective target function locally.

10. The method of claim 9, wherein the derivation of optima in step (6) comprises using kernel density estimate to assess in which areas there are enough simulations to provide sufficient coverage and using kernel regression to estimate objective target function locally in all other areas, wherein the parameterization with the highest value based on the kernel regressions is used as the optima.11 . The method of any one of claims 1 to 10, wherein(1) the number of iterations of steps (4) to (6) is determined by setting a threshold value for improvements of the objective target function by each iteration, wherein if said threshold value is not met, no further iteration is carried out; and / or(2) step (7) comprises assessment of the suggested optima by simulating the suggested optima multiple times.

12. The method of any one of claims 1 to 11 , wherein the method is used for design of a breeding program for plants, in particular crop plants, or animals, in particular livestock.

13. Computer-implemented method for the optimization of a breeding program design, comprising the steps of:(1) definition of the multivariate optimization problem and at least one termination criterion, including defining a search space for the parameters to be optimized and defining the breeding objectives (target function) and constraints for the parameters to be optimized;(2) initialization of the first set of parameterizations by creating initial parameterizations within the search space and meeting the constraints to the parameters to be optimized;(3) evaluation of the first set of parameterizations based on the optimization problem defined in step (1) by simulating the respective breeding program with the selected parameterizations by use of stochastic simulation;(4) selection of a number of parameterizations from the area with the highest average value of the objective target function and, optionally, with the highest value of the objective function from (i) the first set of parameterizations based on the evaluation made in step (3) or (ii) from the second or further set of parameterizations based on the evaluation made in the preceding step (6), wherein said number of selected parameterizations is lower than the total number of parameterizations in the set of parameterizations;(5) generation of a second or further set of parameterizations based on the selected parameterizations of step (4) comprising(5a) creating new parameterizations from the selected set of parameterizations obtained in step (4) by combination of selected parameterizations; and(5b) creating new parameterizations from the selected set of parameterizations obtained in step (4) by modification of selected parameterizations;(6) evaluation of the second or further set of parameterizations generated in step (5) based on the optimization problem defined in step (1) and derivation of optima by simulating the respective breeding program with the selected parameterizations by use of stochastic simulation comprising using kernel density estimate to assess in which areas there are enough simulations to provide sufficient coverage and / or using kernel regression to estimate objective target function locally, wherein steps (4) to (6) are repeated until the at least one termination criterion is met; and(7) assessment of the optima derived in step (6) by simulating the suggested optima multiple times.

14. A data processing system comprising means for carrying out at least steps (3) to (6), preferably steps (2) to (6), more preferably steps (1) to (7) of the method of any one of claims 1 to 12.

15. Computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out at least steps (3) to (6), preferably steps (2) to (6), more preferably steps (1) to (7) of the method of any one of claims 1 to 12.

16. Computer-readable data carrier having stored thereon the computer program of claim 14.

Images