Dorian Leroy

Johann Bourcier

Benoit Combemale

When Scientific Software Meets Software Engineering

Keywords:

The development of scientific software relies on the collaboration of various stakeholders for the scientific computing and software engineering activities. Computer languages have an impact on both activities and related concerns, as well as on the engineering principles required to ensure the development of reliable scientific software. The more general-purpose the language is -with low-level, computing-related, system abstractions-the more flexibility it will provide, but also the more rigorous engineering principles and Validation & Verification (V&V) activities it will require from the language user. In this paper, we investigate the different levels of abstraction, linked to the diverse artifacts of the scientific software development process, a software language can propose, and the V&V facilities associated to the corresponding level of abstraction the language can provide to the user. We aim to raise awareness among scientists, engineers and language providers on their shared responsibility in developing reliable scientific software.

SCIENTIFIC SOFTWARE INTO THE WILD

Scientific computing, also known as computational sciences, consists in the use of advanced computing capabilities to understand and solve complex scientific problems (e.g., biological, physical, and social), but also engineering and humanities problems. It aims to predict the behavior or the outcome of a physical system, being natural or man-made. While initially restricted to the activity of certain experts such as scientists or engineers, scientific computing has been made popular and widely exposed to citizens due to major and damaging societal problems such as climate change and more recently pandemics (e.g., . Most of the political Computer decisions regarding these wicked problems are backed up by scientific findings, all based on scientific computing.

Scientific computing is a cross-cutting field, but its heart and soul lie in the development of mathematical models to understand physical systems through their simulations. Those models can be numerical (e.g., systems of differential equations), non-numerical (e.g., agent-based models) or based on analytics (e.g., machine learning models), and capture the behavior of the modeled system. Numerical models can be further refined as continuous or discrete. Simulations of mathematical models correspond to the execution of the computer programs containing these models, the so-called "simulation codes". In this paper, we refer to the subsuming concept of scientific software, which we define as software dedicated to scientific computing and simulation. The development of these scientific software includes both software engineering and scientific computing concerns. Mathematical models and scientific software are therefore tightly intertwined throughout their life cycles. The tools and methods used for their development (e.g., computer languages) have an impact on the definition of both, as well as on the engineering principles required to ensure the development of reliable scientific software.

Due to the ever-increasing intrinsic complexity of such mathematical models and the pursuit of ever more efficient simulations, the scientific computing and software engineering communities have worked together for decades to explore various approaches to assess the reliability of the proposed models and the correctness of their implementations in scientific software [START_REF] Kanewala | Testing scientific software: A systematic literature review[END_REF]. However, when the time comes to implement a new model -and thus new simulation software -scientists and engineers are faced with decisive choices such as what computer language(s) to use (e.g., MATLAB, Mathematica, Fortran, Python, C++, or even an in-house domain-specific language). This choice has important consequences on the expressiveness at hand to implement the model and the corresponding simulation code, but also in terms of software engineering practices to apply to develop reliable and efficient scientific software. As usual, power comes with responsibility! The more general-purpose the language is -with low-level, computing-related, system abstractions -the more flexibility and performance it may provide, but also the more rigorous engineering principles and Validation & Verification (V&V) activities it will require from the language user to develop a reliable piece of scientific software.

However, most scientists and engineers do not have a background in software engineering and are thus not aware of its best practices beyond programming (e.g., version management, component reuse, continuous integration, unit testing) [START_REF] Wilson | Best practices for scientific computing[END_REF], which has actually prompted initiatives aiming to tackle the issue, such as the Research Software Engineering [START_REF] Cohen | The four pillars of research software engineering[END_REF]. Since the final goal is to build and apply the model encoded in the simulation code, the code itself is merely a means to that end. Final stakeholders (e.g., citizens, policy and decision makers, research institutions, system users...) may even be ignorant of the importance of code for science and engineering. This situation has recently been brought to light by Konrad Hinsen in a post reporting criticisms about an epidemics simulation code (http://blog.khinsen.net/posts/2020/05/18/anopen-letter-to-software-engineers-criticizingneil-ferguson-s-epidemics-simulation-code/), and the subsequent discussions on the blog as well as on social media (https://twitter.com/ khinsen/status/1262307434632282112).

In this paper, we explore the overall scientific software development process, we provide an integrated view of the scientific computing and software engineering activities, artifacts and roles, and we discuss the trade-offs on the required V&V activities with regard to the computer languages at hand. This work offers an holistic view through which we introduce the existing V&V approaches from the literature to support such activities and the accountability of the respective roles. We aim to raise awareness among scientists, engineers and language providers on their shared responsibility in developing reliable scientific software. mon classification separates numerical models (e.g., systems of differential equations) from nonnumerical models (e.g., agent-based models). We focus our discussion on scientific software based solely on continuous numerical models, which we will address as mathematical models in the rest of this paper.

ON THE ENGINEERING OF SCIENTIFIC SOFTWARE

Software engineers generally follow an engineering process, aka. Software Development Life Cycle (SDLC), that introduces the required activities to produce software. For instance, the V-model (cf. green part in Figure 1) has been introduced to highlight the relationships between each activity of the development life cycle and its associated activity of V&V. In the rest of the paper, we use the V-model as a support for the identification of the main software engineering activities, regardless of the associated method (e.g., incremental, iterative or agile). The same activities can be retrieved in all methods, possibly with different time dimensions (e.g., a V-model in each sprint of an agile method).

In Figure 1, we propose an engineering process adapted to scientific software development in the form of a V-model, aka. scientific V-model, subsuming the original software engineering (SE) V-model. We introduce the roles, activities and artifacts of this V-model, with an emphasis on V&V concerns. The V-model is composed of two branches: the first one on the left is top-down, corresponding to the successive elaboration of the various artifacts, and the second on the right is bottom-up, representing the validation process of the scientific software. At each level, a V&V activity on the upward branch is located across from an artifact on the downward branch. The purpose of the V&V activity is to validate the scientific software with respect to the artifact in question. This is represented by the dashed arrows in Figure 1.

The development of scientific software involves several activities that manipulate different types of artifacts. Thus, the design of scientific software based on mathematical models requires the involvement and cooperation of various stakeholders ranging from scientists and engineers to experts in software engineering or numerical analysis. These stakeholders play one of three roles (according to the development context, one person might endorse more than one role): scientists as domain experts (S in Figure 1), numerical analysts as experts in the discretization of a continuous phenomenon on a computer (N in Figure 1), and software engineers as experts of software development to deliver the expected services (SE in Figure 1). It is worth noting the very different background and expertise of these three roles, and the tight collaboration required across the development process.

The tasks linked to each role are associated with the activities and artifacts forming the Vmodel represented in Figure 1. From a set of observations of a phenomenon and theories (illustrated in the upper left part of the figure), scientists develop a mathematical model, a.k.a. continuous numerical model, often taking the form of systems of differential equations. Since these systems of equations cannot generally be solved analytically, they must be solved numerically, which requires a discretization process. This consists in applying a discretization method (e.g., finite differences, finite volumes) to obtain a numerical scheme. This numerical scheme is a discretized mathematical model, which specifies the sequence of computations to execute, given a discretization of the domain according to discretization parameters (e.g., space, time), an initial state (e.g., initial temperature in every point of the domain) and inputs from the environment (e.g., rainfall), to compute a numerical (as opposed to analytical) solution to the model. The implementation of the discretization process requires specific skills and it is undertaken by numerical analysts (N in dark grey in Figure 1), in collaboration with scientists who act as domain experts (S in light grey in Figure 1).

From that point on, software engineering concerns enter the development process, as well as the software engineers. While numerical analysts are required as domain expert (N in light grey in Figure 1), such software engineering concerns require specific skills brought by software engineers (SE in dark grey in Figure 1). This is represented on the middle lower part of Figure 1 with a second, nested V-model (in green) corresponding to the usual V-model used in software engineering. This V-model effectively bridges the gap between the numerical scheme -which can be considered as the stakeholder requirements for scientific software -and the software engineering process. It encompasses all the software engineering concerns to obtain a reliable working software, including performance, concurrency, targeted architectures, memory management, data access, and so on.

On the left top-down branch, the stakeholder requirements are refined and extended with more specific software concerns corresponding to the different levels of organization of the software. The elaboration of the different specifications follows a standard process. First, the software requirements are formulated from the global stakeholder requirements and they specify the expectations for the software as a whole. They contain information about the general software organization and the expected behaviors of the software in various scenarios of interest (e.g., use case diagrams). Then, the general design is drawn to more specifically describe the different components that shape the software (e.g., component diagrams). This leads to the formulation of the detailed design of the software which aims to document all the structural details of the software (e.g., class and sequence diagrams) in order to then allow the implementation of the software in a clear and efficient way. Using the various requirements contained in the different previous specifications, the implementation of the scientific software is generated.

The resulting scientific software is then validated according to the stakeholder requirements, through V&V activities under the form of unit, integration, system, and acceptance testing, all encompassed by implementation testing from the point of view of the numerical analyst. First, unit testing ensures that the implemented scientific software respects the guidelines of the detailed design, meaning that each unit of the software code performs as expected (e.g., mock-based testing). The fulfillment of the general design specifications are checked through integration testing as it makes sure that the different components interact correctly (e.g., Big Bang testing, incremental testing). Then, system testing controls that the implementation meets the software requirements by validating the complete and fully integrated software -considered here as the system of interest -thanks to various testing techniques (e.g., usability testing, load testing, regression testing, functional testing). Finally, compliance with the stakeholders' expectations is validated through acceptance testing (e.g., user acceptance testing, business acceptance testing, operational testing).

Once the software engineering V&V concerns have been addressed, the conformity of the scientific software to its original mathematical model is ensured throughout a step called discretization testing. During this testing activity, numerical analysts aim to validate whether the scientific software produces results that are within an acceptable margin with regards to the results

4

Computer that would be obtained by analytically solving the mathematical model. The techniques used for discretization testing vary in their level of rigor, depending mostly on the availability of an exact solution to the mathematical model. Without such a solution, techniques such as symmetry, conservation, or Galilean invariance testing can be employed [START_REF] Robertson | Computational Eulerian hydrodynamics and Galilean invariance[END_REF].

These techniques leverage domain knowledge that the numerical solution should exhibit specific characteristics such as symmetry when provided with symmetric geometry and conditions, conservation of some physical properties (e.g., mass, energy), or be unaffected by operational changes such as inverting the axes of the coordinate system. Note that, as these tests rely on domain knowledge, some might not be used from one field to another (e.g., when no symmetry is present). Another technique that can be used without an exact solution is code-to-code comparison [START_REF] Oberkampf | On the role of code comparisons in verification and validation[END_REF], which consists in running a simulation of the same mathematical model embedded in another piece of scientific software, and comparing the results.

Yet, obtaining the highest degree of confidence in the discretized model requires techniques leveraging an exact solution to the mathematical model. When an analytical solution is not available, exact solutions can be obtained through the use of the method of manufactured solutions [START_REF] Roache | Code verification by the method of manufactured solutions[END_REF], which constrains the simulator to compute a solution chosen beforehand (the socalled manufactured solution). Whether analytical or manufactured, an exact solution allows to employ rigorous techniques such as discretization error quantification, iterative convergence testing, and order-of-accuracy testing, respectively in order of increasing rigor [START_REF] Oberkampf | Verification and validation in scientific computing[END_REF]. In essence, these techniques aim to assess whether the discretization error, i.e., the difference between the numerical and the exact solutions, tends toward zero as the value of the discretization parameters decreases, and thus that the numerical scheme is a correct discretization of the mathematical model.

Finally, the last V&V step we designate as model testing is conducted by scientists. Model testing (a.k.a. model validation) aims to statistically quantify the disagreement between the numerical solutions and the available experimental data, and to assess whether the resulting uncertainty is acceptable within the application domain of the model [START_REF] Roy | A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing[END_REF]. Ideally, these data are collected through a set of validation experiments. Validation experiments are a special kind of experiment aiming to measure the conditions of the experiment as precisely and completely as possible, with less emphasis on controlling the environment than for other kinds of experiments [START_REF] Oberkampf | Assessment criteria for computational fluid dynamics model validation experiments[END_REF]. This distinction allows to run simulations of the model whose inputs match very closely those of the experiments, and thus to assess the fidelity, in these precise conditions, of the outputs of the simulations with regard to collected measurements.

A general approach to model testing is proposed in [START_REF] Roy | A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing[END_REF]. Broadly, this approach consists in (i) identifying and characterizing mathematically every source of uncertainty (aleatory and epistemic) [START_REF] Helton | Representation of analysis results involving aleatory and epistemic uncertainty[END_REF] in both the scientific software and experimental measurements used as a reference, (ii) propagating those uncertainties to the outputs through techniques such as sensitivity analysis [START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF], (iii) defining a validation metric, i.e., a mathematical operator comparing the numerical solution and the experimental measurements, and applying it to produce an assessment of the validity of the model over the validation domain, and finally (iv) employing statistical techniques such as regression fits of the validity metric over the validation domain, which enables to interpolate or extrapolate the validity metric over the domain of intended use, thereby providing the uncertainty of the simulation results over the domain of application.

Overall, the different testing activities of the scientific V-model can be related to the corresponding steps of the software engineering Vmodel: the implementation testing step can be conceptually compared to the unit and integration testing step, the discretization testing to the system testing, and the model testing to the final scientific software acceptance testing.

These activities serve the same purpose but use different techniques, according to the nature of the artifact under test. The observations fulfill the role of the stakeholders requirements in that they directly represent the essence of the phenomenon to be modeled. The mathematical model can be considered as the software requirements as it encompasses the global scientific behavior of the physical system the software should encode. The numerical scheme reflects the purpose of the global and detailed designs as it specifies how the numerical solution to the mathematical model should be computed in the software, and how each step of computation is related to one another, in a similar way as the behavior and interaction of units and components are described in the design specifications.

However, the techniques used are distinctive as they depend on very different natures of artifacts -present on the left top-down branch of the V-model depicted in Figure 1 -on which they validate the software. For instance, discretization testing, which we relate to system testing in the context of scientific software development, aims at validating the scientific software with regards to the mathematical model. Since the mathematical model consists of mathematical equations, associated V&V techniques manipulate mathematical concepts such as the convergence numerical sequences in the case of iterative convergence testing. Alternatively, the system testing phase of the nested V-model uses the software requirements as a reference to assess the validity of the software. As those are software specifications, the testing techniques manipulate software states and behaviors. They ensure that the different states are consistent from a software perspective. For example, that multiple simultaneous runs of the software do not affect the output of each individual run.

POWER COMES WITH RESPONSIBILITY!

The overall development process is complex and challenging when it comes to ensure reliable scientific software. To mitigate this complexity and allow scientists and numerical analysts to directly implement their scientific software, existing computer languages come with dedicated constructs that balance the responsibility between the language user and the language provider. We explore various computer languages at hand to implement scientific software according to their provided constructs, and discuss the consequences on the expressive power and on the required V&V activities. We illustrate this in Figure 2, where we rely on the V-cycle for the development of scientific software proposed in Figure 1. In Figure 2, the orange and blue colors are used to identify who supplies what artifacts and who is responsible for what V&V activities. The orange color refers to language users, while the blue color refers to language providers.

REUSING EXISTING SIMULATORS

Figure 2a shows the situation where scientists directly use existing scientific software such as MODFLOW [START_REF] Harbaugh | MODFLOW-2005, the US Geological Survey modular ground-water model: the groundwater flow process[END_REF] or MITgcm [START_REF] Adcroft | MITgcm documentation[END_REF]. In this situation, language users are scientists configuring the models encoded into an existing piece of scientific software for their specific case. The language used in this case corresponds to a configuration language used to configure simulators and run simulations according to specific parameters (e.g., FloPy for MODFLOW [START_REF] Bakker | Scripting MODFLOW model development using Python and FloPy[END_REF]). Language users specify fixed input quantities such as the geometry, or physical modeling parameters of the system and obtain a simulator specific to a certain context (e.g., groundwater flows in a specific watershed area). Then, for each simulation, language users specify the input quantities that vary from one simulation to another, such as the initial and boundary conditions, or the system excitation over time.

Oftentimes, some of those input quantities are unknown (e.g., when they cannot be measured in the field), or known only with a high uncertainty. In such cases, language users will need to conduct model calibration to reduce the parametric uncertainty of their model and obtain useful results. Over the years, a number of techniques have been developed to conduct model calibration [START_REF] Jiang | A sequential calibration and validation framework for model uncertainty quantification and reduction[END_REF]. Yet, model calibration cannot eliminate every source of uncertainty. Thus, some means of propagating uncertainties over the inputs through the model to obtain the corresponding output uncertainties are required (e.g., sensitivity analysis [START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF]). While numerous black-box methods exist, white-box (a.k.a. open-box) methods that leverage domain knowledge receive more attention [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: A data-driven, physics-informed Bayesian approach[END_REF]. However, such methods require access to internal variables of the embedded model. The latter must thus be exposed as configurable simulation parameters to enable the use of certain state-of-the-art techniques.

In the situation depicted in Figure 2a, the simulator offers abstractions specific to the embedded mathematical model (e.g., soil perme-6 Figure 2: Scientific software development: responsibilities among the language user (orange) and the language provider (blue), depending on the abstraction level of the provided language constructs. ability, porosity), which users can parameterize to run context-specific simulations through the underlying scientific software. Thus, developers of a simulator must first endorse all the responsibilities associated with the development of a complete piece of scientific software, as covered in the previous section, including implementation, discretization, and model testing.

Computer

Stakeholders Requirements

SE V-Model

An important language-related responsibility that befalls simulator developers is to clearly identify, document, and enforce as much as possible the validity envelope of the abstractions exposed by the configuration language of the simulator, to prevent its misuse. In the general sense, we define the validity envelope of a language as the set of valid programs one can write with that language. In the case of a configuration language dedicated to a specific piece of scientific software (such as FloPy), this validity envelope must prevent unsound module and parameter combinations, and enforce the value of parameters to be within their bounds (e.g., non-negative volume of rainfall).

Due to combinatorial explosion, an abstraction of the domain is generally required to identify and enforce the validity envelope of configuration languages, as far as unsound combinations of parameters are concerned particularly. This will, in turn, necessarily result in false positives or false negatives. Techniques such as sensitivity analysis can be used to identify the relationships between the different parameters, and thus detect such unsound combinations.

LANGUAGES FOR MATHEMATICAL MODELS

When existing scientific software do not embed the desired mathematical model, one will need to define their own. Defining a mathematical model and deriving the corresponding scientific software can be done directly, using languages such as Mathematica (https://www. wolfram.com/mathematica) or MATLAB (https: //matlab.mathworks.com). Figure 2b depicts the development of such languages and the development of scientific software using them.

In this scenario, language users are scientists implementing their mathematical model using a language that offers constructs at the corresponding abstraction level (i.e., continuous mathematics), and allows them to derive the corresponding piece of scientific software. Thus, they endorse the responsibility of addressing the V&V concerns associated with model testing -identified in the previous section -since those are out of scope of such a language. In addition, in order to validate the implementation of their mathematical model with regard to its specification, language users perform discretization testing.

However, as explained in the previous section, when no analytical solution to the mathematical model is available, the most rigorous techniques for discretization testing (e.g., iterative convergence testing, order-of-accuracy testing) can only be used by employing the method of manufactured solutions [START_REF] Roache | Code verification by the method of manufactured solutions[END_REF], which requires to inject source terms into the discretized model. Therefore, to enable the use of these techniques, the language needs to either provide tools allowing it, or to directly give access to the discretized model. Otherwise, only less rigorous techniques such as those described in the previous section can be employed to conduct discretization testing.

Language providers are in charge of designing continuous mathematical constructs (e.g., differential blocks in MATLAB's block diagrams) allowing language users to define their mathematical models, and also developing the infrastructure (e.g., compiler) to derive reliable working scientific software from any mathematical model. This language infrastructure performs automated and possibly configurable discretization and scientific software derivation (as indicated by the blue numerical scheme and scientific software boxes on Figure 2b). So, language providers must handle the corresponding V&V concerns (i.e., the concerns pertaining to the discretization of their language constructs), as well as the software engineering concerns -encapsulated into the SE V-model -which are tied to the implementation of their discretized language constructs.

For that reason, in addition to software engineers, numerical analysts take an active part in language design, and scientists come in assistance to provide the language requirements and criteria for the verification and validation of the language constructs and of their discretization. While the techniques used to address the V&V concerns related to the discretization of language constructs and of plain mathematical models are the same, they are here applied to individual language constructs and combinations thereof, instead of complete mathematical models. The objective is to cover as many valid uses of the language as possible to check that its validity envelope subsumes the application domain intended for the language.

Finally, to allow language users to validate the implementation of their mathematical model, tools for discretization testing should be provided, in particular tools allowing an automated use of the method of manufactured solutions. This allows language users to perform the required validation at the level of abstraction of the language (the continuous mathematical level in this case), which is an essential functionality of software languages.

LANGUAGES FOR NUMERICAL SCHEMES

When more control is required over what discretization methods to apply and how to apply them, and when languages at the mathematical level of abstraction do not fit this need, one has to handle the discretization step oneself. Figure 2c depicts a situation in which numerical analysts use languages dedicated to the definition of numerical schemes (or with the right abstractions to do so) such as NabLab [START_REF] Lelandais | Fostering metamodels and grammars within a dedicated environment for HPC: the NabLab environment (tool demo)[END_REF], Julia [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF], SciPy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF], GNU Octave (https://www.gnu.org/ software/octave/), or Blitz++ (https://github.com/ blitzpp/blitz). These languages allow to automatically derive the corresponding piece of scientific software, on which the scientific V&V activities listed in the previous section can be conducted, without having to handle software engineering concerns.

In this scenario, language users are numerical analysts who apply discretization methods on mathematical models provided by scientists to obtain numerical schemes, that they implement using a language thaht offers constructs at the corresponding level of abstraction, and that allows them to derive a piece of scientific software 8

Computer for the discretized model. Therefore, they endorse the responsibility of addressing the V&V concerns corresponding to discretization testing while scientists endorse the one associated with model testing. Indeed, those responsibilities are out of the scope of a language at the numerical scheme level of abstraction, as indicated by the orange color of the corresponding boxes in Figure 2c. Finally, language users must validate their implementation of their numerical scheme with regard to its specification, which is the aim of implementation testing.

In this situation, implementation testing is performed with the aid of tool support from the language that allows language users to stay at the same abstraction level (i.e., discrete mathematics), without having to consider system concerns such as memory management or concurrency.

Language providers expose discrete mathematical constructs that allow language users to implement their numerical scheme. From the numerical schemes defined with those constructs, the infrastructure of the language (e.g., model transformations, interpreters, compilers, code generators) derives the corresponding piece of scientific software. Therefore, language providers must tackle the system concerns (e.g., memory management, concurrency, data management) and corresponding V&V concerns of the language infrastructure to make it possible to tackle the concerns of the scientific software thus generated. In essence, this means that the concerns must be addressed for any valid numerical scheme written by the language user, which is depicted in Figure 2c by the blue color of the SE V-model and scientific software.

While this task is mainly undertaken by software engineers, the assistance of numerical analysts is required to both define the required expressivity for the language (which in turn allows to derive the stakeholder requirements from a given numerical scheme), and perform the necessary V&V activities to ensure that the effective semantics of the derived piece of software is equivalent to the semantics of the provided numerical schemes. In other words, software engineers require the assistance of numerical analysts to tackle concerns related to the abstractions at the discrete mathematics level, and clearly define the validity envelope of the language.

Then, to foster the robustness of the V&V process spanning the entire V-model, language providers should offer tool support allowing language users to perform implementation testing at the discrete mathematics level of abstraction. This tooling can take the form of a dedicated debugger, and optimally includes facilities dedicated to the analysis of numerical schemes such as plot drawing, monitoring variable values and watch expressions (i.e., exogenous expressions of no interest for actual simulations, evaluated for debugging purposes only), and conditional breakpoints and step-by-step execution to detect and investigate faults in the numerical scheme.

LANGUAGES FOR SCIENTIFIC SOFTWARE

Finally, when performance is critical to a particular piece of scientific software, be it in terms of memory or execution time, specific hardware and architectures need to be targeted (e.g., heterogeneous and performance-oriented infrastructures such as distributed CPU-GPU architectures). If no language for either mathematical models or numerical schemes supports the targeted type of architecture, system-level languages such as C, C++, and Fortran can be used, conjointly with APIs such as OpenMP, and standards such as MPI and SYCL. Continuous progress in the hardware domain and in the field of high performance computing (HPC) pushes the developers of performance-critical scientific software to use such system-level abstractions, as they are often the only way to leverage the latest HPC-specific technologies.

Figure 2d depicts developments where language users express the particularities of their simulator with regard to all the concerns involved in the development of scientific software, ranging from the mathematical model, to the encoded numerical scheme, to system-level concerns such as concurrency, memory, and data handling. In this case, language users must also endorse the responsibilities of addressing the V&V concerns corresponding to each step in the development process, from the V&V of the aforementioned system concerns, to discretization testing, to model testing. In such a situation, another possibility for language users is to take the role of language provider for an existing language for mathematical models or numerical schemes. This places them as language providers in the situations depicted in Figure 2c or Figure 2b. As such, they can broaden the range of architectures that can be targeted by the existing language by including their architecture of choice, and thus foster the use of the language in the community as a language at a higher abstraction level than system-level languages that nonetheless allows to target state-of-the-art hardware.

KEY TAKEAWAYS AND LOOKING AHEAD

When developing scientific software, dedicated languages for mathematical models and numerical schemes (e.g., MATLAB, NabLab, SciPy) bring benefits to users from the scientific community with regard to languages at a low abstraction level (e.g., C++, Fortran), as low-level V&V and performance concerns are addressed as part of the development of those high-level, dedicated languages. This means that users of such high-level languages, i.e., scientists and numerical analysts, can leave aside software engineering concerns to focus on their area of expertise and the associated V&V concerns, resulting in an increased robustness of the end-to-end development process of scientific software.

In exchange for these benefits to language users, designers of high-level, dedicated languages face additional challenges: the higher the abstraction level offered by a language is, the more steps and complexity in the V&V process of the language itself there are. These languages have to guarantee the correctness and performance of the scientific software resulting from any valid piece of code, mathematical model, or numerical scheme written with them, and provide tools to language users to validate their use of the language. For instance, developing a dedicated language at the mathematical level requires the designers to ensure not only that their implementation and discretization testing cover any mathematical model written with the language, but also that users are able to perform discretization testing specific to their precise piece of scientific software.

We highlight in this paper the importance of extensively documenting the V&V processes that are conducted as part of the development of a language. This allows language users to clearly identify what V&V concerns are or not already addressed as part of a language. Consequently, users can pick languages according to both their goals in terms of scientific software development, and to the V&V concerns they are able to continuously address themselves according to the state of the current knowledge. This way, scientists are less likely to have to handle software engineering concerns when they do not want to double as software engineers. In addition, this allows to capitalize on the aforementioned development and V&V overhead that comes with the development of those languages.

Another key point is the need for language providers to offer, as part of the language infrastructure, facilities to allow language users to conduct V&V activities tailored to their specific simulation model. This enables the continuity of the V&V process, from user-defined program to scientific software, and fosters both transparency of the complete development process, and trust in the resulting scientific software.

The work presented in this paper leads to a clarification of the roles of the various stakeholders taking part in the development of scientific software. A reasoned approach for the development of reliable scientific software has been proposed that allows to clearly characterize the validity envelope of this type of software. We believe that the generalized use of our approach will allow to systematically characterize the validity envelope of scientific software, to make it explicit and thus lead to a better and safer use of these software.

ACKNOWLEDGMENT

We would like to thank our colleagues at the Institut de Recherche en Informatique et Systèmes Aléatoires (IRISA), at the Geoscience department of the Observatoire des Sciences de l'Univers de Rennes (OSUR), as well as at the Commissariat à l'énergie atomique et aux énergies alternatives (CEA), for their input and providing such interesting conversations that contributed to the writing and improvement of this paper. We also want to thank Hinsen Konrad and Gordon Blair for their valuable comments, which helped to improve the paper.

10

Computer

Figure 1 :

 1 Figure 1: Overall scientific software development process across the V-model.

  Scientific software development using an existing simulator. (eg., Modflow[START_REF] Harbaugh | MODFLOW-2005, the US Geological Survey modular ground-water model: the groundwater flow process[END_REF]) Scientific software development using languages for mathematical models. (eg., Matlab, Mathematica) Scientific software development using languages for numerical schemes. (eg., Blitz++, SciPy, NabLab) Scientific software development using languages with system abstractions. (eg., C++, Python)

Published by the IEEE Computer Society © 2021 IEEE

Scientific software can rely on various types of scientific models that represent the behavior of the physical system under study. A com-

May/June 2021

Dorian Leroy is a PhD student doing an international PhD between the JKU Linz (Austria) and the University of Rennes 1 (France), currently in the DiverSE team in Rennes. His research interests lie in the field of Software Language Engineering and include metaprogramming approaches and generic V&V facilities. Contact him at dorian.leroy@inria.fr.

June


Sébastien Verel
email: verel@i3s.unice.fr

Gabriela Ochoa

Marco Tomassini
email: marco.tomassini@unil.ch

The Connectivity of NK Landscapes' Basins: A Network Analysis

We propose a network characterization of combinatorial fitness landscapes by adapting the notion of inherent networks proposed for energy surfaces (3). We use the well-known family of N K landscapes as an example. In our case the inherent network is the graph where the vertices represent the local maxima in the landscape, and the edges account for the transition probabilities between their corresponding basins of attraction. We exhaustively extracted such networks on representative small N K landscape instances, and performed a statistical characterization of their properties. We found that most of these network properties can be related to the search difficulty on the underlying N K landscapes with varying values of K.

Introduction

Local optima are the very feature of a landscape that makes it rugged. Therefore, an understanding of the distribution of local optima is of utmost importance for the understanding of a landscape. Combinatorial landscapes refer to the finite search spaces generated by important discrete problems such as the traveling salesman problem and many others. A property of some combinatorial landscapes, which has been often observed, is that on average, local optima are much closer to the optimum than are randomly chosen points, and closer to each other than random points would be. In other words, the local optima are not randomly distributed, rather they tend to be clustered in a "central massif" (or "big valley" if we are minimising). This globally convex landscape structure has been observed in the N K family of landscapes (6), and in other combinatorial optimization problems, such as the traveling salesman problem (2), graph bipartitioning (7), and flowshop scheduling (9).

In this study we seek to provide fundamental new insights into the structural organization of the local optima in N K landscapes, particularly into the connectivity of their basins of attraction. Combinatorial landscapes can be seen as a graph whose vertices are the possible configurations. If two configurations can be transformed into each other by a suitable operator move, then we can trace an edge between them. The resulting graph, with an indication of the fitness at each vertex, is a representation of the given problem fitness landscape. A useful simplification of the graphs for the energy landscapes of atomic clusters was introduced in (3; 4). The idea consists in taking as vertices of the graph not all the possible configurations, but only those that correspond to energy minima. For atomic clusters these are well-known, at least for relatively small assemblages. Two minima are considered connected, and thus an edge is traced between them, if the energy barrier separating them is sufficiently low. In this case there is a transition state, meaning that the system can jump from one minimum to the other by thermal fluctuations going through a saddle point in the energy hyper-surface. The values of these activation energies are mostly known experimentally or can be determined by simulation. In this way, a network can be built which is called the "inherent structure" or "inherent network" in (3). We propose a network characterization of combinatorial fitness landscapes by adapting the notion of inherent networks described above. We use the well-known family of N K landscapes as an example because they are a useful tunable benchmark that can provide interesting information for more realistic combinatorial landscapes. In our case the inherent network is the graph where the vertices are all the local maxima and the edges account for transition probabilities between their corresponding basins of attraction. We exhaustively extract such networks on representative small N K landscape instances, and perform a statistical characterization of their properties. Our analysis was inspired, in particular, by the work of (3; 4) on energy landscapes, and in general, by the field of complex networks (8). The study of networks has exploded across the academic world since the late 90's. Researchers from the mathematical, biological, and social sciences have made substantial progress on some previously intractable problems, bringing new techniques, reformulating old ideas, and uncovering unexpected connections between seemingly different problems. We aim here at bringing the tools of network analysis for the study of problem hardness in combinatorial optimization. The next section describes how combinatorial landscapes are mapped onto networks, and includes the relevant def-initions and algorithms used in our study. The empirical network analysis of our selected N K landscape instances is presented next, followed by our conclusions and ideas for future work.

Landscapes as Networks

To model a physical energy landscape as a network, (4) needed to decide first on a definition both of a state of the system and how two states were connected. The states and their connections will then provide the nodes and edges of the network. For systems with continuous degrees of freedom, the author achieved this through the 'inherent structure' mapping. In this mapping each point in configuration space is associated with the minimum (or 'inherent structure') reached by following a steepest-descent path from that point. This mapping divides the configuration space into basins of attraction surrounding each minimum on the energy landscape.

Our goal is to adapt this idea to the context of combinatorial optimization. In our case, the nodes of the graph can be straightforwardly defined as the local maxima of the landscape. These maxima are obtained exhaustively by running a best-improvement local search algorithm (HillClimbing, see Algorithm 1) from every configuration of the search space. The definition of the edges, however, is a much more delicate matter. In our initial attempt (10) we considered that two maxima i and j were connected (with an undirected and unweighted edge), if there exists at least one pair of solutions at Hamming distance one s i and s j , one in each basin of attraction (b i and b j ). We found empirically on small instances of N K landscapes, that such definition produced densely connected graphs, with very low (≤ 2) average path length between nodes for all K. Therefore, apart from the already known increase in the number of optima with increasing K, no other network property accounted for the increase in search difficulty. Furthermore, a single pair of neighbors between adjacent basins, may not realistically account for actual basin transitions occurring when using common heuristic search algorithms. These considerations, motivated us to search for an alternative definition of the edges connecting local optima. In particular, we decided to associate weights to the edges that account for the transition probabilities between the basins of attraction of the local optima. More details on the relevant algorithms and formal definitions are given below.

Definitions and Algorithms

Definition: Fitness landscape. A landscape is a triplet (S, V, f ) where S is a set of potential solutions i.e. a search space, V : S -→ 2 S , a neighborhood structure, is a function that assigns to every s ∈ S a set of neighbors V (s), and f : S -→ R is a fitness function that can be pictured as the height of the corresponding solutions.

In our study, the search space is composed by binary strings of length N , therefore its size is 2 N . The neighborhood is defined by the minimum possible move on a binary search space, that is, the 1-move or bit-flip operation. In consequence, for any given string s of length N , the neighborhood size is |V (s)| = N .

The HillClimbing algorithm to determine the local optima and therefore define the basins of attraction, is given below:

Algorithm 1 HillClimbing Choose initial solution s ∈ S repeat choose s ′ ∈ V (s) such that f (s ′ ) = max x∈V (s) f (x) if f (s) < f (s ′ ) then s ← s ′ end if until s is a Local optimum Definition: Local optimum. A local optimum is a solution s * such that ∀s ∈ V (s * ), f (s) < f (s * ).
The HillClimbing algorithm defines a mapping from the search space S to the set of locally optimal solutions S * .

Definition: Basin of attraction. The basin of attraction of a local optimum i ∈ S is the set b i = {s ∈ S | HillClimbing(s) = i}. The size of the basin of attraction of a local optima i is the cardinality of b i .

Definition: Edge weight. Notice that for a non-neutral fitness landscapes, as are N K landscapes, the basins of attraction as defined above, produce a partition of the configuration space S. Therefore, S = ∪ i∈S * b i and ∀i ∈ S ∀j = i, b i ∩ b j = ∅ For each solutions s and s ′ , let us define p(s → s ′ ) as the probability to pass from s to s ′ with the bit-flip operator. In the case of binary strings of size N , and the neighborhood defined by the bit-flip operation, there are N neighbors for each solution, therefore:

if s ′ ∈ V (s) , p(s → s ′ ) = 1 N and if s ′ ∈ V (s) , p(s → s ′ ) = 0.
We can now define the probability to pass from a solution s ∈ S to a solution belonging to the basin b j , as:

p(s → b j ) = s ′ ∈bj p(s → s ′ ) Notice that p(s → b j ) ≤ 1.
Thus, the total probability of going from basin b i to basin b j is the average over all s ∈ b i of the transition probabilities to solutions s ′ ∈ b j :

p(b i → b j ) = 1 ♯b i s∈bi p(s → b j )
♯b i is the size of the basin b i . We are now prepared to define our 'inherent' network or network of local optima. Definition: Local optima network. The local optima network G = (S * , E) is the graph where the nodes are the local optima1 , and there is an edge e ij ∈ E with the weight

w ij = p(b i → b j ) between two nodes i and j if p(b i → b j ) > 0.
According to our definition of weights, w ij = p(b i → b j ) may be different than w ji = p(b j → b i ). Two weights are needed in general, and we have an oriented transition graph.

Empirical Basin and Network Analysis

The N K family of landscapes ( 6) is a problem-independent model for constructing multimodal landscapes that can gradually be tuned from smooth to rugged. In the model, N refers to the number of (binary) genes in the genotype (i.e. the string length) and K to the number of genes that influence a particular gene (the epistatic interactions). By increasing the value of K from 0 to N -1, N K landscapes can be tuned from smooth to rugged. The k variables that form the context of the fitness contribution of gene s i can be chosen according to different models. The two most widely studied models are the random neighborhood model, where the k variables are chosen randomly according to a uniform distribution among the n -1 variables other than s i , and the adjacent neighborhood model, in which the k variables that are closest to s i in a total ordering s 1 , s 2 , . . . , s n (using periodic boundaries). No significant differences between the two models were found in (6) in terms of global properties of the respective families of landscapes, such as mean number of local optima or autocorrelation length. Similarly, our preliminary studies on the characteristics of the N K landscape optima networks did not show noticeable differences between the two neighborhood models. Therefore, we conducted our full study on the more general random model.

In order to avoid sampling problems that could bias the results, we used the largest values of N that can still be analysed exhaustively with reasonable computational resources. We thus extracted the local optima networks of landscape instances with N = 14, 16, 18, and K = 2, 4, 6, ..., N -2, N -1. For each pair of N and K values, 30 randomly generated instances were explored. Therefore, the networks statistics reported below represent the average behaviour of 30 independent instances.

Basins of Attraction

Besides the maxima network, it is useful to describe the associated basins of attraction as these play a key role in search algorithms. Furthermore, some characteristics of the basins can be related to the optima network features. The notion of the basin of attraction of a local maximum has been presented before. We have exhaustively computed the size and number of all the basins of attraction for N = 16 and N = 18 and for all even K values plus K = N -1. In this section, we analyze the basins of attraction from several points of view as it is described below.

Global optimum basin size versus K. In Fig. 1 we plot the average size of the basin corresponding to the global maximum for N = 16 and N = 18, and all values of K studied. The trend is clear: the basin shrinks very quickly with increasing K. This confirms that the higher the K value, the more difficult for an stochastic search algorithm to locate the basin of attraction of the global optimum Number of basins of a given size. Fig. 2 shows the cumulative distribution of the number of basins of a given size (with regression line) for a representative instances with N = 18, K = 4 . Table 1 shows the average (of 30 independent landscapes) correlation coefficients and linear regression coefficients (intercept ( ᾱ) and slope ( β)) between the number of nodes and the basin sizes for instances with N = 16, 18. Notice that distribution decays exponentially Table 1: Correlation coefficient ( ρ), and linear regression coefficients (intercept ( ᾱ) and slope ( β)) of the relationship between the basin size of optima and the cumulative number of nodes of a given (basin) size ( in logarithmic scale: log(p(s)) = α + βs + ǫ). The average and standard deviation values over 30 instances, are shown. or faster for the lower K and it is closer to exponential for the higher K. This could be relevant to theoretical studies that estimate the size of attraction basins (see for example (5)). These studies often assume that the basin sizes are uniformly distributed, which is not the case for the N K landscapes studied here. From the slopes β of the regression lines (table 1) one can see that high values of K give rise to steeper distributions (higher β values). This indicates that there are fewer basins of large size for large values of K. Basins are thus broader for low values of K, which is consistent with the fact that those landscapes are smoother.

N = 16 K ρ ᾱ β 2 -0.

Fitness of local optima versus their basin sizes.

The scatter-plots in Fig. 3 illustrate the correlation between the basin sizes of local maxima (in logarithmic scale) and their fitness values. Two representative instances for N = 18 and K = 4, 8 are shown. Notice that there is a clear positive correlation between the fitness values of maxima and their basins' sizes. In other words, the higher the peak the wider tend to be its basin of attraction. Therefore, on average, with a stochastic local search algorithm, the global optimum would be easier to find than any other local optimum. This may seem surprising. But, we have to keep in mind that as the number of local optima increases (with increasing K), the global optimum basin is more difficult to reach by a stochastic local search algorithm (see Fig. 1). This observation offers a mental picture of N K landscapes: we can consider the landscape as composed of a large number of mountains (each corresponding to a basin of attraction), and those mountains are wider the taller the hilltops. Moreover, the size of a mountain basin grows exponentially with its hight.

General Network Statistics

We now briefly describe the statistical measures used for our analysis of maxima networks.

The standard clustering coefficient (8) does not consider weighted edges. We thus use the weighted clustering measure proposed by (1), which combines the topological information with the weight distribution of the network:

c w (i) = 1 s i (k i -1) j,h w ij + w ih 2 a ij a jh a hi where s i = j =i w ij , a nm = 1 if w nm > 0, a nm = 0 if w nm = 0 and k i = j =i a ij .
For each triple formed in the neighborhood of the vertex i, c w (i) counts the weight of the two participating edges of the vertex i. C w is defined as the weighted clustering coefficient averaged over all vertices of the network.

The standard topological characterization of networks is obtained by the analysis of the probability distribution p(k) that a randomly chosen vertex has degree k. For our weighted networks, a characterization of weights is obtained by the connectivity and weight distribution p(w) that any given edge has weight w.

In our study, for each node i, the sum of weights from the node i is equal to 1. So, an important measure is the weight w ii of self-connecting edges (remaining in the same node). We have the relation: w ii + s i = 1. The vertex strength, s i , is defined as s i = j∈V(i)-{i} w ij , where the sum is over the set V(i) -{i} of neighbors of i (1). The strength of a node is a generalization of the node's connectivity giving information about the number and importance of the edges.

Another network measure we report here is disparity (1) Y 2 (i), which measures how heterogeneous is the contributions of the edges of node i to the total weight (strength):

Y 2 (i) = j =i w ij s i 2
The disparity could be averaged over the node with the same degree k. If all weights are nearby of s i /k, the disparity for nodes of degree k is nearby 1/k.

Finally, in order to compute the average distance (shortest path) between two nodes on the optima network of a given landscape, we considered the expected number of bit-flip mutations to pass from one basin to the other. This expected number can be computed by considering the inverse of the transition probabilities between basins. In other words, if we attach to the edges the inverse of the transition probabilities, this value would represent the average number of random mutations to pass from one basin to the other. More formally, the distance (expected number of bit-flip mutations) between two nodes is defined by d ij = 1/w ij where w ij = p(b i → b j ). Now, we can define the length of a path between two nodes as being the sum of these distances along the edges that connect the respective basins.

Detailed Study of Network Features

In this section we study in more depth some network features which can be related to stochastic local search difficulty on the underlying fitness landscapes. Table 2 reports the average (over 30 independent instances for each N and K) of the network properties described. nv and ne are, respectively, the mean number of vertices and the mean number of edges of the graph for a given K rounded to the next integer. Cw is the mean weighted clustering coefficient. Ȳ is the mean disparity, and d is the mean path length.

Clustering Coefficients. The fourth column of table 2 lists the average values of the weighted clustering coefficients for all N and K. It is apparent that the clustering coefficients decrease regularly with increasing K for all N .

For the standard unweighed clustering, this would mean that the larger K is, the less likely that two maxima which are connected to a third one are themselves connected. Taking weights, i.e. transition probabilities into account this means that either there are fewer transitions between neighboring basins for high K, and/or the transitions are less likely to occur. This confirms from a network point of view the common knowledge that search difficulty increases with K. Shortest Path to the Global Optimum. The average shortest path lengths d are listed in the sixth column of table 2. Fig. 4 (top) is a graphical illustration of the average shortest path length between optima for all the studied N K landscapes. Notice that the shortest path increases with N , this is to be expected since the number of optima increases exponentially with N . More interestingly, for a given N the shortest path increases with K, up to K = 10, and then it stagnates and even decreases slightly for the N = 18. This is consistent with the well known fact that the search difficulty in N K landscapes increases with K. However, some paths are more relevant from the point of view of a stochastic local search algorithm following a trajectory over the maxima network. In order to better illustrate the relationship of this network property with the search difficulty by heuristic local search algorithms, Fig. 4 (bottom) shows the shortest path length to the global optimum from all the other optima in the landscape. The trend is clear, the path lengths to the optimum increase steadily with increasing K.

Weight Distribution Here we report on the weight distributions p(w) of the maxima network edges. Fig. 5 shows the empirical probability distribution function for the cases N = 16 and N = 18 (logarithmic binning has been used on the x-axis). The case N = 14 is similar but is not reported here because it is much more noisy for K = 2 and 4 due to the small size of the graphs in these cases (see table 2). One can see that the weights, i.e. the transition probabilities between neighboring basins are small. The distributions are far from uniform and, for both N = 16 and N = 18, the low K have longer tails. For high K the decay is faster. This seems to indicate that, on average, the transition probabilities are higher for low K. Disparity Fig. 6 depicts the disparity coefficient as defined in the previous section for N = 16, 18. An interesting observation is that the disparity (i.e. dishomogeneity) in the weights of a node's outcoming links tends to decrease steadily with increasing K. This reflects that for high K the transitions to other basins tend to become equally likely, which is another indication that the landscape, and thus its representative maxima network, becomes more random and difficult to search.

When K increases, the number of edges increases and the number of edges with a weight over a certain threshold increases too (see fig. 5). Therefore, for small K, each node is connected with a small number of nodes each with a relative high weight. On the other hand, for large K, the weights be-come more homogeneous in the neighbourhood, that is, for each node, all the neighboring basins are at similar distance.

If we suppose that edges with higher weights are likely to be connected to nodes with larger basins (an intuition that we need to confirm in future work). Then, as the larger basins tend to have higher fitness (see Fig. 3), the path to higher fitness values would be easier to find for lower K than for larger K. Boundary of basins. Fig. 7 shows the averages, over all the nodes in the network, of the weights w ii (i.e the probabilities of remaining in the same basin after a bit-flip mutation). Notice that the weights w ii are much higher when compared to those w ij with j = i (see Fig. 5). In particular, for K = 2, 50% of the random bit-flip mutations will produce a solution within the same basin of attraction. These average probabilities of remaining within the same basin, are above 12% for the higher values of K. Notice that the averages are nearly the same regardless the value of N , but decrease with the epistatic parameter K.

The exploration of new basins with the random bit-flip mutation seems to be, therefore, easier for large K than for low K. But, as the number of basins increases, and the fitness correlation between neighboring solutions decreases with increasing K, it becomes harder to find the global maxima for large K. This result suggests that the dynamic of stochastic local search algorithms on N K landscapes with large K is different than that with lower values of K, with the former engaging in more random exploration of basins.

The boundary of a basin of attraction can be defined as the set of configurations within a basin that have at least one neighbor's solution in another basin. Conversely, the interior of a basin is composed by the configurations that have all their neighbors in the same basin. Table 3 gives the average number of configurations in the interior of basins (this statistic is computed on 30 independent landscapes). Notice that the size of the basins' interior is below 1% (except for N = 14, K = 2). Surprisingly, the size of the basins' boundaries is nearly the same as the size of the basins themselves. Therefore, the probability of having a neighboring solution in the same basin is high, but nearly all the solutions have a neighbor solution in another basin. Thus, the interior basins seem to be "hollow", a picture which is far from the smooth standard representation of landscapes in 2D with real variables where the basins of attraction are visualized as real mountains.

Conclusions

We have proposed a new characterization of combinatorial fitness landscapes using the family of N K landscapes as an example. We have used an extension of the concept of in- herent networks proposed for energy surfaces (3) in order to abstract and simplify the landscape description. In our case the inherent network is the graph where the nodes are all the local maxima and the edges accounts for transition probabilities (using the bit-flip operator) between the local maxima basins of attraction. We have exhaustively obtained these graphs for N = {14, 16, 18}, and for all even values of K, plus K = N -1, and conducted a network analysis on them. Our guiding motivation has been to relate the statistical properties of these networks, to the search diffi- The clustering coefficients suggest that, for high values of K, the transition between a given pair of neighboring basins is less likely to occur. The shortest paths increase with N and, for a given N , they clearly increase with higher K.

The weight distributions indicate that, on average, the transition probabilities are higher for low K.

The disparity coefficients reflect that for high K the transitions to other basins tend to become equally likely, which is an indication of the randomness of the landscape. The construction of the maxima networks requires the determination of the basins of attraction of the corresponding landscapes. We have thus also described the nature of the basins, and found that the size of the basin corresponding to the global maximum becomes smaller with increasing K.

The distribution of the basin sizes is approximately exponential for all N and K, but the basin sizes are larger for low K, another indirect indication of the increasing randomness and difficulty of the landscapes when K becomes large. Furthermore, there is a strong positive correlation between the basin size of maxima and their degrees. Finally, we found that the size of the basins boundaries is roughly the same as the size of basins themselves. Therefore, nearly all the configurations in a given basin have a neighbor solution in another basin. This observation suggests a different landscape picture than the smooth standard representation of 2D landscapes where the basins of attraction are visualized as hilltops.

This study is our first attempt towards a topological and statistical characterization of combinatorial landscapes, from the point of view of complex networks analysis. Much remains to be done. The results should be confirmed for larger instances of N K landscapes. This will require good sampling techniques, or theoretical studies since exhaustive sampling becomes quickly impractical. Other landscape types should also be examined, such as those containing neutrality, which are very common in real-world applications. Finally, the landscape statistical characterization is only a step towards implementing good methods for searching it. We thus hope that our results will help in designing or estimating efficient search techniques and operators.

Figure 1 :Figure 2 :

 12 Figure 1: Average of the relative size of the basin corresponding to the global maximum for each K over 30 landscapes.

Figure 3 :

 3 Figure 3: Correlation between the fitness of local optima and their corresponding basin sizes, for two representative instances with N = 18, K = 4 (top) and K = 8 (bottom).

Figure 4 :

 4 Figure 4: Average distance (shortest path) between nodes (top), and average path length to the optimum from all the other basins (bottom).

Figure 5 :

 5 Figure 5: Probability distribution of the network weights w ij with j = i in logscale on x-axis. Averages of 30 instances for each N and K are reported.

Figure 6 :Figure 7 :

 67 Figure6: Average disparity, Y 2 , of nodes with a given degree k. Average of 30 independent instances for each N and K are reported. The curve 1/k is also reported to compare to random case.

Table 2 :

 2 N K landscapes network properties. Values are averages over 30 random instances, standard deviations are shown as subscripts. n v and n e represent the number of vertexes and edges (rounded to the next integer), Cw , the mean weighted clustering coefficient. Ȳ represent the mean disparity coefficient, d the mean path length (see text for definitions).

	K	nv	ne	Cw	Ȳ	d
			N = 14
	2	146	200131	0.980.0153 0.3670.0934 76194
	4	7010	3163766	0.920.0139 0.1480.0101	896
	6	18415	123271238	0.790.0149 0.0930.0031	1193
	8	35022	258281801	0.660.0153 0.0700.0020	1332
	10	58522	416861488	0.540.0091 0.0580.0010	1391
	12	89622	574201012	0.460.0048 0.0520.0006	1401
	13	108520	65287955	0.420.0045 0.0500.0006	1391
			N = 16
	2	3315	516358	0.960.0245 0.3260.0579	5614
	4	17833	91292930	0.920.0171 0.1370.0111	1268
	6	46029	417914690	0.790.0154 0.0840.0028	1703
	8	89033	933844394	0.650.0102 0.0620.0011	1942
	10	1, 47034	1621394592	0.530.0070 0.0500.0006	2061
	12	2, 25432	2279122670	0.440.0031 0.0430.0003	2071
	14	3, 26429	2907322056	0.380.0022 0.0400.0003	2031
	15	3, 86833	3212032061	0.350.0022 0.0390.0004	2001
			N = 18
	2	5025	15791854	0.950.0291 0.3070.0630	7315
	4	33072	262667056	0.920.0137 0.1270.0081	1749
	6	99473	14644118685	0.780.0155 0.0760.0044	2375
	8	2, 09370	35400918722	0.640.0097 0.0560.0012	2732
	10	3, 61961	62052120318	0.520.0071 0.0440.0007	2921
	12	5, 65759	89974214011	0.430.0037 0.0380.0003	2971
	14	8, 35260	116364011935 0.360.0023 0.0340.0002	2931
	16 11, 79763	14068706622	0.320.0012 0.0320.0001	2831
	17 13, 79577	15247304818	0.300.0009 0.0320.0001	2771

Table 3 :

 3 Average (on 30 independent landscapes for each N and K) of the mean sizes of the basins interiors. .00045 0.0022 0.00012 0.0011 0.00006 10 0.0055 0.00041 0.0031 0.00015 0.0018 0.00006 12 0.0061 0.00041 0.0040 0.00014 0.0025 0.00007 13 0.0059 0.00029 14 0.0045 0.00012 0.0031 0.00004 15 0.0044 0.00014 16 0.0035 0.00005 17 0.0034 0.00006 culty of the underlying combinatorial landscapes when using stochastic local search algorithms (based on the bit-flip operator) to optimize them. We have found clear indications of such relationships, in particular:

	K N = 14	N = 16	N = 18
	2	0.0167 0.02478 0.0050 0.00798 0.0028 0.00435
	4	0.0025 0.00065 0.0012 0.00025 0.0006 0.00010
	6	0.0029 0.00037 0.0014 0.00015 0.0007 0.00009
	8	0.0043 0

Since each maximum has its associated basin, G also describes the interconnection of basins.

Acknowledgements. Thanks to Professor Edmund Burke for useful advice and support. This work was partially funded by EPSRC grant number EP/D061571/1.


Amélie Barozet

Pablo Chacón

Juan Cortés
email: juan.cortes@laas.fr

Current approaches to flexible loop modeling

Keywords: Loop modeling, Protein flexibility, Conformational sampling, Structure prediction, Energy landscapes

Loops are key components of protein structures, involved in many biological functions. Due to their conformational variability, the structural investigation of loops is a difficult topic, requiring a combination of experimental and computational methods. This paper provides a brief overview of current computational approaches to flexible loop modeling, and presents the main ingredients of the most standard protocols. Despite great progress in recent years, accurately modeling the conformational variability of long flexible loops remains a challenging problem. Future advances in this field will likely come from a tight coupling of experimental and computational techniques, which would enable a better understanding of the relationships between loop sequence, structural flexibility, and functional roles. In fine, accurate loop modeling will open the road to loop design problems of interest for applications in biomedicine and biotechnology.

such as molecular recognition, allosteric regulation or signaling. On average, loops correspond to about one third of a protein sequence, and up to one half of the residues of enzyme active sites are located in loops [START_REF] Bartlett | Analysis of catalytic residues in enzyme active sites[END_REF].

Protein loops can have very different dynamic behaviors. In some cases, loops adopt a relatively stable conformation in the folded state of the protein, being part of the overall protein scaffold. Based on this static perspective, some efforts have been made for the classification of loop structures [START_REF] Oliva | An automated classification of the structure of protein loops[END_REF]. However, conformational changes in loop regions are frequently observed to drive enzyme function or regulation [START_REF] Malabanan | A role for flexible loops in enzyme catalysis[END_REF]. For example, flexible loops (P-loop and T-loop) modulate protein kinases activity. Loops are also essential in many protein-ligand/cofactor interactions. For instance, extracellular loops of G-Protein-Coupled Receptors (GPCRs) are crucial for ligand recognition and binding [START_REF] Hilger | Structure and dynamics of GPCR signaling complexes[END_REF]. Note that GPCRs and protein kinases are by far the major drug targets for pharmaceutical industry.

Together with their variable motion amplitude, loop conformational transitions can take place in a broad range of timescales. A comprehensive list of cases demonstrating the functional importance of protein loops can be found in a recent review [START_REF] Papaleo | The role of protein loops and linkers in conformational dynamics and allostery[END_REF]. Figure 1 is an illustrative summary of loop functions and of the interplay between their structural properties and functional roles.

To better understand the dynamic nature of loops, a suitable representation of their conformational space is based on the concept of energy landscape [START_REF] Wales | Energy landscapes: Applications to clusters, biomolecules and glasses[END_REF]. Although one loop can exist in multiple conformations, they are not equally likely to be adopted. More precisely, the probability of each possible state depends on its associated free energy. The energy landscape provides an overall view by displaying the energy at each point of the conformational space. The topology of this landscape gives insight into the main conformational and dynamic properties of the loop: the basins constitute meta-stable conformations, while saddle regions are associated with transition states. A visual representation can then be obtained by projecting this landscape along two meaningful coordinates, but other representations based on stationary points of the landscape can also be utilized. This is illustrated in Figure 2.

The structural investigation and characterization of flexible loops is challenging. X-ray crystallography, which is the most widely used experimental method to determine high-resolution protein structures, only provides static snapshots in particular experimental conditions (e.g. presence/absence of a ligand). Furthermore, flexible regions are often missing in X-ray structures, or even in electron microscopy maps, because of lacking electron density. According to [START_REF] Djinovic-Carugo | Missing strings of residues in protein crystal structures[END_REF], 69% of the structures deposited in the Protein Data Bank (PDB) have missing fragments, and this percentage increases up to 80% for structures solved with a resolution above 2.0 Å. In more than 90% of the cases, these missing fragments are located in loops or unstructured terminal regions. Other experimental techniques, such as X-ray solution scattering [START_REF] Petoukhov | Addition of missing loops and domains to protein models by x-ray solution scattering[END_REF] or nuclear magnetic resonance (NMR) [START_REF] Boehr | The dynamic energy landscape of dihydrofolate reductase catalysis[END_REF], can provide very relevant but limited structural and kinetic information regarding these flexible regions. . Therefore, since an accurate atomistic representation of the diverse conformations adopted by the loop from experimental measurements is extremely difficult, computational methods are an essential complement to study flexible protein loops.

From a practical point of view, loop modeling techniques usually implement a multi-stage protocol, which varies depending on the underlying modeling approach. Figure 3 aims to provide a general scheme in which most of the existing loop modeling techniques can fit. In brief, modeling can be divided into three major stages: (1) conformational sampling or search, (2) scoring and clustering, and (3) a final post-processing or refinement stage.

The core of loop modeling is the exploration of the vast conformational space to be performed at the sampling/search step. Note that the volume of the conformational space grows exponentially with loop length. Very different approaches to protein loop modeling have been proposed over the years depending on the sampling/search protocol [START_REF] Shehu | Modeling structures and motions of loops in protein molecules[END_REF][START_REF] Papaleo | The role of protein loops and linkers in conformational dynamics and allostery[END_REF][START_REF] Kundert | Computational design of structured loops for new protein functions[END_REF]. Overall, they can be divided into three main classes: knowledge-based, ab initio and hybrid approaches. Knowledge-based approaches, also called template-based or homology-based, exploit structural repositories to retrieve observed loop conformations for a given set of sequence and geometric descriptors about the anchoring points [START_REF] Choi | FREAD revisited: Accurate loop structure prediction using a database search algorithm[END_REF][START_REF] Messih | LoopIng: a Template-Based Tool for Predicting the Structure of Protein Loops[END_REF][START_REF] Shirvanizadeh | Loop modelling 1.0[END_REF][START_REF] Karami | DaReUS-Loop: Accurate loop modeling using fragments from remote or unrelated proteins[END_REF]. These methods are in general computationally fast since they do not rely on expensive simulations, but they are limited by the availability of suitable loop conformations from known protein structures, and only provide a small set of solutions.

Especially for long loops, there is not enough structural data currently to cover the whole conformational space.

On the other end of the spectrum, ab initio approaches can perform a wider sampling of the conformational space, for example by exhaustively varying the torsional angles of the loop [START_REF] Jacobson | A hierarchical approach to all-atom protein loop prediction[END_REF][START_REF] López-Blanco | RCD+: Fast loop modeling server[END_REF]. Their computational cost is highly variable, mainly depending on the underlying sampling method, but overall they are more demanding than knowledge-based approaches. To a lesser extent, they are also limited when dealing with long loops, as sampling coverage decreases with loop length.

Numerous loop modeling methods combine ideas from knowledge-based and ab initio approaches. Hybridization is performed at different levels: some methods apply a consensus between results obtained from different types of approaches [START_REF] Deane | CODA: A combined algorithm for predicting the structurally variable regions of protein models[END_REF], while other methods apply a more intricate combination of strategies to explore the conformational space. In particular, many hybrid approaches use small fragments from structural databases within an ab initio sampling technique [START_REF] Stein | Improvements to robotics-inspired conformational sampling in Rosetta[END_REF][START_REF] Tang | Fast protein loop sampling and structure prediction using distance-guided sequential chain-growth Monte Carlo method[END_REF][START_REF] Park | Protein loop modeling using a new hybrid energy function and its application to modeling in inaccurate structural environments[END_REF][START_REF] Marks | Sphinx: Merging knowledge-based and ab initio approaches to improve protein loop prediction[END_REF][START_REF] Barozet | A reinforcement-learning-based approach to enhance exhaustive protein loop sampling[END_REF]. The interest of one type of approach over the others may depend on the specific application and on the nature of the loop. For instance, knowledge-based approaches can be a suitable choice to predict the most probable conformations of more structurally conserved L1, L2, L3, H1, and H2 antibody CDR loops, whereas the more exhaustive sampling performed by ab initio and hybrid approaches will be more adequate to the highly variable H3 loop.

One of the main ingredients of ab initio and hybrid sampling approaches are loop-closure methods. They are required to ensure that N-terminal and Cterminal ends of the loop keep bonded to the anchoring points in the protein, while correct bond geometry is preserved. They are usually based on the concept of inverse kinematics, very common in robotics. Two classes of loop-closure methods have been proposed: those based on numerical optimization techniques [START_REF] Canutescu | Cyclic coordinate descent: A robotics algorithm for protein loop closure[END_REF][START_REF] Chys | Random coordinate descent with spinor-matrices and geometric filters for efficient loop closure[END_REF], and those based on analytical solutions [START_REF] Coutsias | A kinematic view of loop closure[END_REF][START_REF] Cortés | Geometric algorithms for the conformational analysis of long protein loops[END_REF]. Numerical methods are conceptually simpler and directly applicable to long loops, but analytical methods are successfully used within iterative algorithms for the same purpose.

Scoring methods are another key ingredient, required to evaluate the quality of the sampled loop states. Scoring methods can also be classified into knowledge-based (López-Blanco and Chacón, 2019; [START_REF] Dong | Optimized atomic statistical potentials: assessment of protein interfaces and loops[END_REF], physicsbased (generally based on molecular mechanics approaches, not quantum physics models) and hybrid methods [START_REF] Park | Protein loop modeling using a new hybrid energy function and its application to modeling in inaccurate structural environments[END_REF]Alford et al., 2017). As before, knowledge-based scoring functions, typically linked to a coarse-grain representation of the loop, are the most computationally efficient alternative. However, their implicit simplifications prevent a faithful representation of the rugged energy surface of the all-atom conformational space. Conversely, the computational cost of detailed physics-based energy function prevents a thorough exploration of the space. Hybrid approaches try to balance the best of these two worlds. A common practice is to use less accurate scoring functions in an initial search, and an atomistic detailed potential to refine only the high-scoring loop states. A suitable combination of sampling protocols and scoring functions is the major burden for loop modeling, and the key for future progress in the field.

In addition to the aforementioned loop modeling methods, Monte Carlo (MC) and molecular dynamics (MD) simulations can be used to investigate thermodynamic and kinetic properties of loops. Note however that in this context, MC-derived approaches are mainly employed as a sampling engine within more sophisticated loop prediction techniques, e.g. [START_REF] Stein | Improvements to robotics-inspired conformational sampling in Rosetta[END_REF][START_REF] Tang | Fast protein loop sampling and structure prediction using distance-guided sequential chain-growth Monte Carlo method[END_REF]. Enhanced MD methods [START_REF] Bernardi | Enhanced sampling techniques in molecular dynamics simulations of biological systems[END_REF] are suitable tools to study loop dynamics [START_REF] Papaleo | The role of protein loops and linkers in conformational dynamics and allostery[END_REF]. For example, MD simulations coupled with a higher-level representation such as Markov state models are a relevant approach to investigate the role of loop conformational changes in enzyme catalysis [START_REF] Liao | Loop motion in triosephosphate isomerase is not a simple open and shut case[END_REF]. Nevertheless, due to their very high computational cost, MD simulations are often only applied in a post-processing stage to refine solutions provided by much faster loop sampling methods [START_REF] Lee | Effective protein model structure refinement by loop modeling and overall relaxation[END_REF].

Figure 4 provides a (non-exhaustive) list of loop modeling protocols, related techniques for sampling and scoring, and structural repositories. Note also that other methods (in addition to MD) have been proposed to model loop conformational transitions rather that loop conformations, but they are not covered by this review.

Despite huge progress over the last two decades, flexible loop modeling remains an open problem [START_REF] Marks | Predicting loop conformational ensembles[END_REF][START_REF] Barozet | Protein loops with multiple meta-stable conformations: a challenge for sampling and scoring methods[END_REF]. State-of-theart methods are able to predict stable conformations of relatively short loops (up to 12 residues). However, accurately sampling, scoring and representing the great diversity of loop conformations and transitions between them is still a computational challenge, in particular for long loops. Nowadays, Machine/Deep-Learning (ML/DL) approaches are showing huge potential in many areas, including structural bioinformatics [START_REF] Gao | Deep learning in protein structural modeling and design[END_REF][START_REF] Gkeka | Machine learning force fields and coarse-grained variables in molecular dynamics: Application to materials and biological systems[END_REF][START_REF] Pakhrin | Deep learning-based advances in protein structure prediction[END_REF]. In the context of loop modeling, ML/DL-based approaches have already provided good results as scoring methods [START_REF] Ruffolo | Geometric potentials from deep learning improve prediction of CDR H3 loop structures[END_REF]. To the best of our knowledge, the ability of ML/DL methods to generate meaningful conformational ensemble models of flexible loops remains to be demonstrated, but it is undoubtedly a promising research direction.

The main limitation for the development of more accurate and general loop modeling methods is the lack of experimental data. As mentioned before, loop flexibility is a challenge for biophysical methods: X-ray crystallography provides only static and possibly biased snapshots, NMR methods have difficulties dealing with large proteins, other methods (such as small-angle X-ray scattering or Förster resonance energy transfer) only provide coarse-grained constraints to build models. Therefore, new integrated approaches, tightly coupling several experimental and computational methods, are necessary for advances in this field.

Future improvements in loop modeling will help us to better understand the relationships between loop sequence, structural flexibility and biological function. Meanwhile, current approaches, with their limitations, are ready to use in challenging problems such as loop design in enzymes and antibodies [START_REF] Kundert | Computational design of structured loops for new protein functions[END_REF]. Encouraging results in these fields bring us closer to the ultimate goal of designing loops able to recognize partner molecules with high specificity and affinity, or to enhance the catalytic activity of enzymes thanks to optimized dynamic properties.

Loop Modeling Protein Loops

The importance of protein loop modeling

Conformational variability

Loop Structural Flexibility In some cases, they fluctuate around a main conformational state. For example, some CDR loops in antibodies adopt a single conformation suited to the recognition of a specific antigen. In other cases, a loop can adopt two main conformational states. For instance, this may correspond to a loop opening and closing the access to the active site of an enzyme. In general, the energy landscape of a loop is more complex, involving several basins corresponding to different meta-stable states. These multiple states can also be functionally important, for example in the case of multi-specific antibodies.

Figure 1 :

 1 Figure1: Illustration of important functional roles of protein loops. Their location at the protein surface together with their conformational versatility make them key actors in intermolecular interactions. Their flexibility constitutes a challenge for their structural investigation, which requires a combination of experimental and computational methods.

Figure 2 :

 2 Figure 2: Loops exhibit different levels of conformational variability.In some cases, they fluctuate around a main conformational state. For example, some CDR loops in antibodies adopt a single conformation suited to the recognition of a specific antigen. In other cases, a loop can adopt two main conformational states. For instance, this may correspond to a loop opening and closing the access to the active site of an enzyme. In general, the

Acknowledgments

This work has been supported by the French National Research Agency (ANR) through grant ANR-19-PI3A-0004 and by the Spanish grants BFU2016-76220-P and PID2019-109041GB-C21 (AEI,10.13039/501100011033).

Authorship contribution statement

All the authors contributed equally to this work.

Declaration of competing interest

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

Alford, R., Leaver-Fay, A., Jeliazkov, J.R., O'Meara, M.J., DiMaio, F.P., Park, H., Shapovalov, M.V., Renfrew, P.D., Mulligan, V.K., Kappel, K., Labonte, J.W., Pacella, M.S., Bonneau, R., Bradley, P., Dunbrack, R., Das, R., Baker,

Conformational Sampling/Search Methods

Target Loop

Sequence