est destinée au dépôt

Gambierol and analogs reverse paralysis caused by d-tubocurarine and botulinum type-A neurotoxin at vertebrate neuromuscular junctions

Jordi Molgó, Sébastien Schlumberger, Makoto Sasaki, Haruhiko Fuwa, M.

Carmen Louzao, Luis M Botana, Denis Servent, Evelyne Benoit

To cite this version:

Jordi Abstract: Gambierol is a marine polycyclic ether toxin, first isolated with ciguatoxins by Satake, Murata, and Yasumoto (1993) from cultured Gambierdiscus toxicus dinoflagellates, collected in French Polynesia. The chemical synthesis of gambierol and truncated analogues has been achieved successfully allowing detailed analyses of their mode of action. Interestingly, gambierol and analogs inhibit voltage-gated K + (Kv) channels. It is well known that Kv channels play important roles in modulating acetylcholine (ACh) release from motor nerve terminals of vertebrate neuromuscular junctions. Therefore, the aim of our study was to investigate the action of synthetic gambierol and analogues on K + currents in nerve cells and nerve terminals, using electrophysiological techniques. In addition, gambierol was studied on neuromuscular junctions in which muscle nicotinic ACh receptors have been blocked with d-tubocurarine (postsynaptic block), or in junctions in which quantal ACh release has been greatly reduced (presynaptic block) by botulinum neurotoxin type-A (BoNT/A). This neurotoxin induces a selective proteolytic cleavage of SNAP-25 an important nerve terminal protein for the neurotransmitter release process. Our results showed that gambierol and analogues inhibited Kv channels in neuronal cells in a concentration-dependent manner. Also, nanomolar concentrations of gambierol prolonged the duration of the presynaptic action potential in motor nerve terminals, as revealed by focal current recording. This indicates that Kv channels in nerve terminals were also blocked by the polyether. The loading of the calcium-sensitive dye Fluo-3/AM to frog motor nerve terminals and the use of spectrometric imaging techniques showed that gambierol increased the transient calcium fluorescence signals in response to 1-10 Hz nerve-stimulation, indicating that blockade of Kv channels enhanced calcium entry into nerve terminals. At the neuromuscular junction, such actions of gambierol are essential for increasing nerve-evoked quantal ACh release and can explain the observed reversal of the neuromuscular block produced by d-tubocurarine, a competitive inhibitor of the muscle nicotinic ACh receptor. Interestingly, gambierol also increased the quantal content of endplate potentials in BoNT/A-intoxicated mouse neuromuscular junctions. These results indicate that blockade of Kv channels and enhanced calcium entry into terminals by gambierol increased the probability of quantal ACh release, supporting the view that SNAP-25 cleaved by BoNT/A can enter into the protein-complex required for calciumdependent neurotransmitter release. In conclusion, our results suggest that gambierol and analogues can have potential medical application in neuromuscular pathologies in which it is necessary to antagonize pre-or post-synaptic neuromuscular blockade, or both.

Disclosure of Interest: None Declared

  Molgó, Sébastien Schlumberger, Makoto Sasaki, Haruhiko Fuwa, M. Carmen Louzao, et al.. Gambierol and analogs reverse paralysis caused by d-tubocurarine and botulinum type-A neurotoxin at vertebrate neuromuscular junctions. ICHA 2018 -18th International conference on Harmful algae, Institut Francais de Recherche pour l'Exploitation de la Mer (Ifremer), Oct 2018, Nantes, France. reverse paralysis caused by d-tubocurarine and botulinum type-A neurotoxin at vertebrate neuromuscular junctions J. Molgó 1 2,* , S. Schlumberger 2 , M. Sasaki 3 , H. Fuwa 3 , M. C. Louzao 4 , L. M. Botana 4 , D. Servent 1 , E. Benoit 1 2 1 Institut des sciences du vivant Frédéric Joliot, Service d'Ingénierie Moléculaire des Protéines, Université Paris-Saclay, CEA, 2 Institut des Neurosciences Paris-Saclay, UMR 9197 CNRS / Université Paris-Sud, CNRS, Gif sur Yvette, France, 3 Graduate School of Life Sciences, Tohoku University, Sendai, Japan, 4 Departamento de Farmacología, Facultad de Veterinaria, Universidad de Santiago de Compostela, Lugo, Spain


Steeven Janny
email: steeven.janny@insa-lyon.fr

Madiha Nadri
email: madiha.nadri-wolf@univ-lyon1.fr

Julie Digne
email: julie.digne@cnrs.fr

SPACE AND TIME CONTINUOUS PHYSICS SIMULATION FROM PARTIAL OBSERVATIONS

Modern techniques for physical simulations rely on numerical schemes and meshrefinement methods to address trade-offs between precision and complexity, but these handcrafted solutions are tedious and require high computational power. Data-driven methods based on large-scale machine learning promise high adaptivity by integrating long-range dependencies more directly and efficiently. In this work, we focus on fluid dynamics and address the shortcomings of a large part of the literature, which are based on fixed support for computations and predictions in the form of regular or irregular grids. We propose a novel setup to perform predictions in a continuous spatial and temporal domain while being trained on sparse observations. We formulate the task as a double observation problem and propose a solution with two interlinked dynamical systems defined on, respectively, the sparse positions and the continuous domain, which allows to forecast and interpolate a solution from the initial condition. Our practical implementation involves recurrent GNNs and a spatio-temporal attention observer capable of interpolating the solution at arbitrary locations. Our model not only generalizes to new initial conditions (as standard auto-regressive models do) but also performs evaluation at arbitrary space and time locations. We evaluate on three standard datasets in fluid dynamics and compare to strong baselines, which are outperformed both in classical settings and in the extended new task requiring continuous predictions.

INTRODUCTION

The Lavoisier conservation principle states that changes in physical quantities in closed regions must be attributed to either input, output, or source terms. By applying this rule at an infinitesimal scale, we retrieve partial differential equations (PDEs) governing the evolution of a large majority of physics scenarios. Consequently, the development of efficient solvers is crucial in various domains involving physical phenomena. While conventional methods (e.g. finite difference or finite volume methods) showed early success in many situations, numerical schemes suffer from high computational complexity, in particular for growing requirements on fidelity and precision. Therefore, there is a need for faster and more versatile simulation tools that are reliable and efficient, and data-driven methods offer a promising opportunity.

Large-scale machine learning offers a natural solution to this problem. In this paper, we address data-driven solvers for physics, but with additional requirements on the behavior of the simulator: R1. Data-driven -the underlying physics equation is assumed to be completely unknown. This includes the PDE, but also the boundary conditions. The dynamics must be discovered from a finite dataset of trajectories, i.e. a collection of observed behaviors from the physical system, R2. Generalization -the method must be capable of handling new initial conditions that do not explicitly belong to the training set, without re-training or fine-tuning, R3. Time and space continuous -the domain of the predicted solution must be continuous in space and time1 so that it can be queried at any arbitrary location within the domain of definition.

These requirements are common in the field but rarely addressed altogether. R1 allows for handling complex phenomena where the exact equation might be unknown, and R2 supports the growing need for faster simulators, which consequently must handle new ICs. Space and time continuity (R3) are also useful properties for standard simulations since the solution can be made as fine as needed in certain complex areas. This task requires learning from sparsely distributed observations only, and without any prior knowledge on the PDE form. In these settings, a standard approach consists of approximating the behavior of a discrete solver, enabling forecasting in an auto-regressive fashion [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF]; [START_REF] Janny | Eagle: Large-scale learning of turbulent fluid dynamics with mesh transformers[END_REF]; [START_REF] Sanchez-Gonzalez | Learning to simulate complex physics with graph networks[END_REF], losing therefore spatial and temporal continuity. Indeed, auto-regressive models assume strong regularities in the data, such as a static spatial lattice and uniform time steps. For these reasons, generalization to new spatial locations or intermediate time steps is not straightforward. These methods satisfy R1 and R2, but not R3. In another trend, Physics-Informed Neural Networks (PINNs) learn a solution on a continuous domain. They leverage the PDE operator to optimize the weights of a neural network representing the solution, and cannot generalize to new ICs, thus violating R1 and R2.

In this paper, we address R1, R2 and R3 altogether in a new setup involving two joint dynamical systems. R1 and R2 are satisfied using an auto-regressive discrete-time dynamics learned from the sparse observations and producing a trajectory in latent space. Then, R3 is achieved with a state observer derived from a second dynamical system in continuous time. This state observer relies on transformer-based cross-attention to enable evaluation at arbitrary spatio-temporal locations. In a nutshell: (a) We propose a new setup to address continuous space and time simulations of physical systems from sparse observation, leveraging insights from control theory. (b) We provide strong theoretical results indicating that our setup is well-suited to address this task compared to existing baselines, which are confirmed experimentally on challenging benchmarks. (c) We provide experimental evidence that our state observer is more powerful than handcrafted interpolations for the targeted task. (d) With experiments on three challenging standard datasets (Navier Yin et al. (2022); [START_REF] Gabriel | On the Effect of the Internal Friction of Fluids on the Motion of Pendulums[END_REF]), Shallow Water Yin et al. (2022); [START_REF] Galewsky | An initial-value problem for testing numerical models of the global shallow-water equations[END_REF], Eagle [START_REF] Janny | Eagle: Large-scale learning of turbulent fluid dynamics with mesh transformers[END_REF], and against state-of-the-art methods (MeshGraphNet (MGN) [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF], DINO Yin et al. (2022), MAgNet [START_REF] Boussif | Magnet: Mesh agnostic neural pde solver[END_REF])), we show that our results generalize to a wider class of problems, with excellent performances.

RELATED WORKS

Autoregressive models -have been extensively used to replicate the behavior of iterative solvers in discrete time, especially in cases where the PDE is unknown or generalization to new initial conditions is needed. These models come in various internal architectures, including convolution-based models for systems observed on a dense uniform grid [START_REF] Stachenfeld | Learned simulators for turbulence[END_REF][START_REF] Le | Disentangling physical dynamics from unknown factors for unsupervised video prediction[END_REF][START_REF] De Bézenac | Deep learning for physical processes: Incorporating prior scientific knowledge[END_REF] and graph neural networks [START_REF] Battaglia | Interaction networks for learning about objects, relations and physics[END_REF] that can adapt to arbitrary spatial discretizations [START_REF] Sanchez-Gonzalez | Learning to simulate complex physics with graph networks[END_REF]Janny et al., 2022a;[START_REF] Li | Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids[END_REF]. Such models have demonstrated a remarkable capacity to produce highly accurate predictions and generalize over long prediction horizons, making them particularly suitable for addressing complex problems such as fluid simulation [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF][START_REF] Han | Predicting physics in meshreduced space with temporal attention[END_REF][START_REF] Janny | Eagle: Large-scale learning of turbulent fluid dynamics with mesh transformers[END_REF]. However, auto-regressive models are inherently limited to a fixed and constant spatio-temporal discretization grid, hindering their capability to evaluate the solution anywhere and at any time. Neural ordinary differential equations (Neural ODE Chen et al. (2018); [START_REF] Dupont | Augmented neural odes[END_REF]) offer a countermeasure to the fixed timestep constraint by learning continuous ODEs on discrete data using an explicit solver, such as Euler or Runge-Kutta methods. In theory, this enables the solution to be evaluated at any temporal location but in practice still relies on the discretization of the time variable. Moreover, extending this approach to PDEs is not straightforward. Contrarily to these approaches, we leverage the auto-regressive capacity and accuracy while allowing arbitrary evaluation of the solution at any point in both time and space.

Continuous solutions for PDEs -date back to the early days of deep learning [START_REF] Mwmg Dissanayake | Neural-network-based approximations for solving partial differential equations[END_REF][START_REF] Isaac E Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF][START_REF] Dimitris | A hybrid neural network-first principles approach to process modeling[END_REF] and have recently experienced a resurgence of interest [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF]2017). Physics-informed neural networks represent the solution directly as a neural network and train the model to minimize a residual loss derived from the PDE. They are mesh-free, which alleviates the need for complex adaptive mesh refinement techniques (mandatory in finite volume methods), and have been successfully applied to a broad range of physical problems [START_REF] Lu | Physics-informed neural networks with hard constraints for inverse design[END_REF][START_REF] George S Misyris | Physics-informed neural networks for power systems[END_REF][START_REF] Zoboli | Deep learning-based output tracking via regulation and contraction theory[END_REF][START_REF] Kissas | Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks[END_REF][START_REF] Yang | Predictive large-eddy-simulation wall modeling via physics-informed neural networks[END_REF][START_REF] Cai | Physicsinformed neural networks (pinns) for fluid mechanics: A review[END_REF], with a growing community proposing architecture designs specifically tailored for PDEs [START_REF] Vincent Sitzmann | Implicit neural representations with periodic activation functions[END_REF][START_REF] Fathony | Multiplicative filter networks[END_REF] as well as new training methods [START_REF] Zeng | Competitive physics informed networks[END_REF][START_REF] Marc | A stable and scalable method for solving initial value pdes with neural networks[END_REF][START_REF] De | Simple initialization and parametrization of sinusoidal networks via their kernel bandwidth[END_REF]). Yet, these models are also known to be difficult to train efficiently [START_REF] Krishnapriyan | Characterizing possible failure modes in physics-informed neural networks[END_REF][START_REF] Wang | When and why pinns fail to train: A neural tangent kernel perspective[END_REF]. Recently, neural operators have attempted to learn a mapping between function space, leveraging kernels in Fourier space (Li et al., 2020b) (FNO) or graphs (Li et al., 2020a) (GNO) to learn the correspondence from the initial condition to the solution at a fixed horizon. While some operator learning frameworks can theoretically generalize to unseen initial conditions and arbitrary locations, we must consider the practical limitations of existing baselines. For instance, FNO requires a static cartesian grid and cannot be directly evaluated outside the training grid. Similarly, GNO can handle arbitrary meshes in theory, but still has limitations in evaluating points outside the training grid and [START_REF] Li | Markov neural operators for learning chaotic systems[END_REF] variant can only be queried at fixed time increments. DeepONet [START_REF] Lu | Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators[END_REF]) can handle free sampling in time and space but is also constrained to a static observation grid.

Continuous and generalizable solvers -represent a significant challenge. Few models satisfy all these conditions. MP-PDE [START_REF] Brandstetter | Message passing neural PDE solvers[END_REF] can handle free-form grids but cannot generalize to different resolutions between train and test, and performs auto-regressive temporal forecasting. Closer to our work, MAgNet [START_REF] Boussif | Magnet: Mesh agnostic neural pde solver[END_REF] proposes to interpolate the observation graph in latent space to new query points before forecasting the solution using graph neural networks. However, they assume prior knowledge of the evaluation mesh and the new query points, use nearest neighbor interpolation instead of trained attention and struggle to generalize to finer grids during test time. In [START_REF] Hua | Efficient continuous spatio-temporal simulation with graph spline networks[END_REF], the auto-regressive MeshGraphNet [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF] is combined with Orthogonal Spline Collocation to allow for arbitrary spatial queries. Finally, DINo [START_REF] Yuan Yin | Continuous pde dynamics forecasting with implicit neural representations[END_REF] proposes a mesh-free, space-time continuous model to address PDE solving. The model uses context adaptation techniques to dynamically adapt the output of an implicit neural representation forward in time. DINo assumes the existence of a latent ODE modeling the temporal evolution of the context vector and learns it as a Neural ODE. In contrast, our method differs from DINo as our model is based on physics forecasting in an auto-regressive manner. We achieve space and time continuity through a learned dynamical attention transformer capable of handling arbitrary locations and points in time. Our design choices allow for generalization on new spatial and temporal locations, ie. not limited to discrete time steps, and new initial conditions while being trainable from sparse observations2 .

CONTINUOUS SOLUTIONS FROM SPARSE OBSERVATIONS

Consider a dynamical system following a Partial Differential Equation (PDE) defined for all (x, t) ∈ Ω × 0, T , with T a positive constant:

ṡ(x, t) = f (s, x, t) ∀(x, t) ∈ Ω × 0, T , s(x, 0) = s 0 (x) ∀x ∈ Ω, s(x, t) = s(x, t) ∀(x, t) ∈ ∂Ω × 0, T (1)
where the state lies in an invariant set s ∈ S, f : S → S is an unknown operator, s 0 : Ω → R n is the initial condition (IC) and s : ∂Ω × 0, T → R n the boundary condition. In what follows, we consider trajectories with shared boundary conditions, hence we omit s from the notation for readability. In practice, the operator f is unknown, and we assume access to a set D of K discrete trajectories from different ICs, s k 0 , sampled at sparse and scattered locations in time and space. Formally, we introduce two finite sets X ⊂ Ω of fixed positions and fixed regularly sampled times T at sampling rate ∆ * . Let S(s 0 , x, t) be the solution of this PDE from IC s 0 , the dataset D is given as: D := S(s k 0 , X , T ) k ∈ 1, K . Our task is formulated as:

Given D, a new initial condition s 0 ∈ S, and a query (x, t) ∈ Ω × 0, T , find the solution of equation 1 at the queried location and from the given IC, that is S(s 0 , x, t). Note that this task involves generalization to new ICs, as well as estimation to unseen spatial locations within Ω and unseen time instants within 0, T . We do not explicitly require extrapolation to instants t > T , although it comes as a side benefit of our approach up to some extent.

THE DOUBLE OBSERVATION PROBLEM

The task implies extracting regularities from weakly informative physical variables that are sparsely measured in space and time, since X and T contain very few elements. Consequently, the possibility to forecast their trajectories from off-the-shelf auto-regressive methods is very unlikely (as confirmed experimentally). To tackle this challenge, we propose an approach accounting for the fact that the phenomenon is not directly observable from the sparse trajectories, but can be deduced from a richer latent state-space in which the dynamics is markovian. We introduce two linked dynamical models lifting sparse observations to dense trajectories guided by observability considerations, namely System 1:

z d [n+1] = f 1 z d [n] s d [n] = h 1 z d [n] , System 2: ṡ(x, t) = f 2 s, x, t z(x, t) = h 2 s, x, t ∀(x, t)∈Ω× 0, T (2)
where for all n ∈ N, we note s d [n] = s(X , n∆) the sparse observation at some instant n∆ (the sampling rate ∆ is not necessarily equal to the sampling rate ∆ * used for data acquisition, which we will exploit during training to improve generalization. This will be detailed later).

System 1 -is a discrete-time dynamical system where the available measurements s d [n] are considered as partial observations of a latent state variable z d [n]. We aim to derive an output predictor from System 1 to forecast trajectories of sparse observations auto-regressively from the sparse IC. As mentioned earlier, sparse observations are unlikely to be sufficient to perform predictions, hence we introduce a richer latent state variable z d in which the dynamics is truly markovian, and observations s d [n] are seen as measurements of the state z d using the function h 1 .

System 2 -is a continuous-time dynamical system describing the evolution of the to-be-predicted dense trajectory S(s 0 , x, t). It introduces continuous observations z(x, t) such that z(X , n∆) = z d [n]. The insight is that the state representation z d [n] obtained from System 1 is designed to contain sufficient information to predict s d [n], but not necessarily to predict the dense state. Formally, z d represents solely the observable part of the state, in the sense of control theory.

At inference time, we forecast at query location (x, t) with a 2-step algorithm: (Step-1) System 1 is used as an output predictor from the sparse IC s d [0], and computes a sequence z[0], z[1], ..., which we refer to as "anchor states". This sequence allows the dynamics to be Markovian, provides sufficient information for the second state estimation step and holds information to predict the sparse observations, allowing supervision during training. (Step-2) We derive a state observer from System 2 leveraging the anchor states over the whole time domain to estimate the dense solution at an arbitrary location in space and time (see figure 1). Importantly, for a given IC, the anchor states are computed only once and reused within System 2 to estimate the solution at different points.

THEORETICAL ANALYSIS

In this section, we introduce theoretical results supporting the use of Systems 1 and 2. In particular, we show that using System 1 to forecast the sparse observations in latent space z d rather than directly operating in the physical space leads to smaller upper bounds on the prediction error. Then, we show the existence of a state estimator from System 2 and compute an upper bound on the estimation error depending on the length of the sequence of anchor states.

Step 1 -consists of computing the sequence of anchor states guided by an output prediction task of the sparse observations. As classically done, we introduce an encoder (formally, a state observer) e s d [0] =z d [0] coupled to System 1 to project the sparse IC into a latent space z d . Following System 1, we compute the anchor states z d auto-regressively (with f 1 ) in the latent space. The sparse observations are extracted from z d using h 1 . In comparison, existing baselines [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF][START_REF] Sanchez-Gonzalez | Learning to simulate complex physics with graph networks[END_REF][START_REF] Stachenfeld | Learned simulators for turbulence[END_REF] maintain the state in the physical space and discard the intermediate latent representation between iterations. Formally, let us consider approximations f1 , ĥ1 , ê (in practice realized as deep networks trained from data D) of f 1 , h 1 and e and compare the prediction algorithm for the classic auto-regressive (AR) approach and ours Classic AR:

ŝar d [n] := ( ĥ1 • f1 • ê) n s d [0] Ours: ŝd [n] := ĥ1 • f1 n • ê s d [0] (3)
Classical AR approaches re-project the latent state into the physical space at each step and repeat "encode-process-decode". Our method encodes the sparse IC, advances the system in the latent space, and decodes toward the physical space at the end. A similar approach has also been explored in [START_REF] Wu | Learning to accelerate partial differential equations via latent global evolution[END_REF]; [START_REF] Kochkov | Learning latent field dynamics of pdes[END_REF], albeit in different contexts, without theoretical analysis.

Proposition 1 Consider a dynamical system of the form of System 1 and assume the existence of a state observer e along with approximations f1 , ĥ1 , ê with Lipschitz constants L f , L h and L e respectively such that

L h L f L e ̸ = 1. If there exist δ f , δ h , δ e ∈ R + such that ∀(z, s) ∈ R nz × R ns |f 1 (z) -f1 (z)| ⩽ δ f , |h 1 (z) -ĥ1 (z)| ⩽ δ h , |e(s) -ê(s)| ⩽ δ e (4)
for the Euclidean norm | • |, then for all integer n > 0, with ŝd [n] and ŝar d [n] as in equation 3,

|s d [n] -ŝd [n]| ⩽ δ h + L h δ f L n f -1 L f -1 + L n f δ e (5) |s d [n] -ŝar d [n]| ⩽ δ L n -1 L -1 (6) with δ = δ h + L h δ f + L h L f δ e and L = L h L f L e . Proof : See appendix B.
This result shows that falling back to the physical space at each time step degrades the upper bound of the prediction error. Indeed, if L < 1, the upper bound converges trivially to zero when n increases, and hence can be ignored. Otherwise, the upper bound for the classic AR scheme appears to be more sensitive to approximation errors δ h , δ f and δ e compared to our approach (for a formal comparison, see appendix C). Intuitively it means that information is lost in the observation space, which thus needs to be re-estimated at each iteration when using the classic AR scheme. By maintaining a state variable in the latent space, we allow this information to flow readily between each step of the simulator (see blue frame in figure 1).

Step 2 -The state estimator builds upon System 2 and relies on the set of anchor states from the previous step to estimate the dense physical state at arbitrary locations in space and time. Formally, we look for a function ψ q leveraging the sequence of anchor states

z d [0], • • • z d [q] (simulated from the sparse IC s d [0]
) to retrieve the dense solution 3 . In what follows, we show that (1) such a function ψ q exists and (2) we compute an upper bound on the estimation error depending on the length of the sequence. To do so, consider the functional which outputs the anchor states from any IC

s 0 ∈ S O p (s 0 )= h 2 s 0 (X ) h 2 S(s 0 , X , ∆) • • • h 2 S(s 0 , X , p∆) = z d [0] • • • z d [p] (7)
In practice, the ground truths z d [n] are not perfectly known, as they are obtained from a data-driven output predictor (step 1) using the sparse IC. Inspired from Janny et al. (2022b), we state:

Proposition 2 Consider a dynamical system defined by System 2 and equation 7. Assume that A1. f 2 is Lipschitz with constant L s , A2. there exists p > 0 and a strictly increasing function α such that ∀s a , s b ∈ S 2 and ∀q

⩾ p O q (s a ) -O q (s b ) ⩾ α(q)|s a -s b | S (8)
where • S is an appropriate norm for S.

Then, ∀q ⩾ p, there exists ψ q such that, for (x, t) ∈ Ω× 0, T and δ n such that ẑd

[n] = z d [n]+δ n , for all n ⩽ q, ψ q z d [0], • • • , z d [q], x, t = S(s 0 , x, t) (9)
S(s 0 , x, t) -ψ q ẑd [0], • • • , z d [q], x, t S ⩽ 2α(q) -1 δ 0|q e Lst . ( 10
)
where

δ 0|q = δ 0 • • • δ q .
Proof: See appendix D. Assumption A2. states that the longer we observe two trajectories from different ICs, the easier it will be to distinguish them, ruling out systems collapsing to the same state. Such systems are uncommon since forecasting their trajectory becomes trivial after some time.

This assumption is related to finite-horizon observability in control theory, a property of dynamical systems guaranteeing that the (markovian) state can be retrieved given a finite number p of past observations. Equation 8 is associated with injectivity of O q , hence the existence of a left inverse mapping the sequence of anchor states to the IC s 0 .

Proposition 2 highlights a trade-off on the performance of ψ q . On one hand, longer sequences of anchor states are harder to predict, leading to a larger |δ 0|q |, which impacts the state estimator ψ q negatively. On the other hand, longer sequences hold more information that can still be leveraged by ψ q to improve its estimation, represented by α(q) -1 in equation 10. In contrast to competing baselines or conventional interpolation algorithms, our approach takes this trade-off into account, by explicitly leveraging the sequence to estimate the dense solution, as will be discussed below.

Discussion and related work -the competing baselines can be analyzed using our setup, yet in a weaker configuration. For instance, one can see Step 2 as an interpolation process, and replace it with a conventional interpolation algorithm, which typically relies on spatial neighbors only. Our method not only exploits spatial neighborhoods but also leverages temporal data, improving the performance, as shown in proposition 2 and empirically corroborated in Section 4.

MAgNet [START_REF] Boussif | Magnet: Mesh agnostic neural pde solver[END_REF] uses a reversed interpolate-forecast scheme compared to ours. The IC s d [0] is interpolated right from the start to estimate s 0 (corresponding to our Step 2, with q=1), and then simulated with an auto-regressive model in the physical space (with the classic AR scheme). Propositions 1 and 2 show that the upper bounds on the estimation and prediction error are higher than ours. Moreover, if the number of query points exceeds the number of known points (|Ω|≫|X |), the input of the auto-regressive solver is filled with noisy interpolations, which impacts performance.

DINo [START_REF] Yuan Yin | Continuous pde dynamics forecasting with implicit neural representations[END_REF]) is a very different approach leveraging a spatial implicit neural representation modulated by a context vector, whose dynamics is modeled via a learned ODE. This approach is radically different than ours and arguably involves stronger hypotheses, such as the existence of a learnable ODE modeling the dynamics of a suitable weight modulation vector. In contrast, our method relies on arguably more sound assumptions, i.e. the existence of an observable discrete dynamics explaining the sparse observation, and the finite-time observability of System 2.

IMPLEMENTATION

The implementation follows the algorithm described in the previous section: (Step-1) rolls out predictions of anchor states from the IC, (Step-2) estimates the state at the query position from these anchor states. We model f1 as a multi-layer Graph Neural Network (GNN) [START_REF] Battaglia | Interaction networks for learning about objects, relations and physics[END_REF]. The anchor states z d [n] are defined at fixed time steps n∆, which might not match ∆ * used in the data T . We found it beneficial to choose ∆=k×∆ * with k>1∈N such that the model can be queried during training on time points t ∈ T that do not match exactly with every time-steps in z d [0], z d [1], ..., but rather on a subset of them, hence encouraging generalization to unseen time. The observation function ĥ1 is an MLP applied on the vector at node level in the graph z d .

The state estimator ψ q is decomposed into a Transformer model [START_REF] Vaswani | Attention is all you need[END_REF] coupled to a recurrent neural network to provide an estimate at query spatio-temporal query position (x, t). First, through cross-attention we translate the set of anchor states z d [n] (one embedding per graph node i and per instant n) into a set of estimates of the continuous variable z(x, t) conditioned at the instant n∆, which we denote z n∆ (x, t) (one embedding per instant n). Following advances in geometric mappings in computer vision [START_REF] Saha | Translating images into maps[END_REF], we use multi-head cross-attention to query from coordinates (x, t) to Keys corresponding to the nodes i in each graph anchor state z d [n], ∀n:

z n∆ (x, t) = f mha Q=ζ ω (x, t), K=V={z d [n] i } + ζ ω (X , n∆) , // attention over nodes i (11)
where Q, K, V are, respectively, Query, Key and Value inputs to the cross-attention layer f mha [START_REF] Vaswani | Attention is all you need[END_REF] and ζ ω a Fourier positional encoding with a learned frequency parameter ω. Finally, we leverage a state observer to estimate the dense solution at the query point from the sequence of conditioned anchor variables, over time. This is achieved with a Gated Recurrent Unit (GRU) [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF] maintaining a hidden state u

[n], u[n] = r gru u[n-1], z n∆ (x, t) , Ŝ(s 0 , x, t) = D (u[q]) , (12)
which shares similarities with conventional state-observer designs in control theory [START_REF] Bernard | Observer design for continuous-time dynamical systems[END_REF]. Finally, an MLP D maps the final GRU hidden state to the desired output, that is, the value of the solution at the desired spatio-temporal coordinate (x, t). See appendix E for details.

TRAINING

Generalization to new input locations during training is promoted by creating artificial generalization situations using sub-sampling techniques of the sparse sets X and T .

Artificial generalization -The anchor states z d [n] are computed at time rate ∆ larger than the available rate ∆ * . This creates situations during training where the state estimator ψ q does not have access to a latent state perfectly matching with the queried time. We propose a similar trick to promote spatial generalization. At each iteration, we sub-sample the (already sparse) IC s d [0] randomly to obtain sd [0] defined on a subset of X . We then compute the anchor states zd using System 1. On the other hand, the query points are selected in the larger set X . Consequently, System 2 is exposed to positions that do not always match with the ones in z d [n]. Note that the complete domain of definition Ω × 0, T remains unseen during training.

Training objective -To reduce training time, we randomly sample M query points (x m , τ m ) in X × T at each iteration, with a probability proportional to the previous error of the model at this point since its last selection (see appendix E) and we minimize the loss

L = K k=1 Lcontinuous M m=1 S(s k 0 , x m , τ m ) -ψ q zd [0|q], x, τ m 2 + Ldynamics ⌊T /∆⌋ n=0 sd [n] -ĥ1 zd [n] 2 , ( 13
) with zd [n] = f1 n • ê sd [0]
. L continuous supervises the model end-to-end, and L dynamics trains the latent anchor states z d to predict the sparse observations from the IC.

EXPERIMENTAL RESULTS

Experimental setup -X ×T results from sub-sampling Ω× 0, T with different rates to control the difficulty of the task. We evaluate on three highly challenging datasets (details in appendix F): Navier [START_REF] Yuan Yin | Continuous pde dynamics forecasting with implicit neural representations[END_REF][START_REF] Gabriel | On the Effect of the Internal Friction of Fluids on the Motion of Pendulums[END_REF] simulates the vorticity of a viscous, incompressible flow driven by a sinusoidal force acting on a square domain with periodic boundary conditions. Shallow Water [START_REF] Yuan Yin | Continuous pde dynamics forecasting with implicit neural representations[END_REF][START_REF] Galewsky | An initial-value problem for testing numerical models of the global shallow-water equations[END_REF] studies the velocity of shallow waters evolving on the tangent surface of a 3D sphere. Eagle [START_REF] Janny | Eagle: Large-scale learning of turbulent fluid dynamics with mesh transformers[END_REF] is a challenging dataset of turbulent airflow generated by a moving drone in a 2D environment with many different scene geometries.

We evaluate our model against three baselines representing the state-of-the-art in continuous simulations. Interpolated MeshGraphNet (MGN) [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF] -3 ).

used auto-regressively and extended to spatiotemporal continuity using physics-agnostic interpolation. MAgNet [START_REF] Boussif | Magnet: Mesh agnostic neural pde solver[END_REF] interpolates the IC at the query position in latent space before using MGN. The original implementation assumes knowledge of the target graph during training, including new queries. When used for superresolution, the authors kept the ratio between the amount of new query points and available points constant. Hence, while MAgNet is queried at unseen locations, it also benefits from more information. In our setup, the model is exposed to a fixed number of points but does not receive more samples during evaluation. This makes our problem more challenging than the one addressed in [START_REF] Boussif | Magnet: Mesh agnostic neural pde solver[END_REF]. DINo (Yin et al., 2022) models the solution as an Implicit Neural Representation (INR) s(x, α t ) where the spatial coordinates x are fed to a MFN [START_REF] Fathony | Multiplicative filter networks[END_REF]) and α t is a context vector modulating the weights of the INR. The dynamics of α is modeled with a Neural-ODE, where the dynamics is an MLP. We share common objectives with DINo and take inspiration from their evaluation tasks yet in a more challenging setup. Details of the baselines are in appendix F. We highlight a caveat on MAgNet: the model can handle a limited amount of new queries, roughly equal to the number of observed points. Our task requires the solution at up to 20 times more queries than available points. In this situation, the graph in MaGNet is dominated by noisy states from interpolation, and the auto-regressive forecaster performs poorly.

During evaluation, we found it beneficial to split the queries into chunks of 10 nodes and to apply the model several times. This strongly improves the performance at the cost of an increased runtime.

Space Continuity -Table 1 compares the spatial interpolation power of our method versus several baselines. The MSE values computed on the training domain (In-X =X ) and outside (Ext-X =Ω\X ) show that our method offers the best performance, especially for the Ext-domain task, which is our aim. To ablate dynamics and evaluate the impact of trained interpolations, we also report the predictions of a Time Oracle which uses sparse ground truth values at all time steps and interpolates (bicubic) spatially. This allows us to assess whether the method is doing better than a simple axiomatic interpolation. While MGN offers competitive in-domain predictions, the cubic interpolation fails to extrapolate reliably on unseen points. This can be seen in the In/Ext gap for Interpolated MGN which is very close to the Time Oracle error. MaGNet, which builds on a similar framework, is hindered by the larger amount of unobserved data in the input mesh. At test time, the same number of initial condition points are provided but the method interpolates substantially more points. DINo achieves a very low In/Ext gap, yet fails on highly (5%) down-sampled tasks. One of the key difference with DINo is that the dynamics relies on an internal ODE for the temporal evolution of a modulation vector. In contrast, our model uses an explicit auto-regressive backbone, and time forecasting is handled in an arguably more meaningful space, which we conjecture to be the reason why we achieve better results (see fig. 5 in the appendix).

Time Continuity -is a step forward in difficulty, as the model needs to interpolate not only to unseen spatial locations (datasets are undersampled at 25%) but also on intermediate timesteps (Ext-T , Table 2). All models perform well on Shallow Water, which is relatively easy. Both DINo and MAgNet leverage a discrete integration scheme (Euler for MAgNet and RK4 for DINo) allowing querying the model between timesteps seen at training. These schemes struggle to capture the data dependencies effectively and therefore the methods fail on Navier (see also Figure 6 for qualitative 1 and2. Table 2: Time Continuity -we evaluate the time interpolation power of our method vs. the baselines. Models are trained and evaluated with 25% of Ω, and with different temporal resolutions (full, half, and quarter of the original). The Spatial Oracle (not comparable!) uses the exact solution at every point in space, and performs temporal interpolation. Evaluation is conducted over 20 frames in the future (10 for Eagle) and we report MSE compared to the ground truth solution (×10 -3 ).

results). Eagle is particularly challenging, the main source of error being the spatial interpolation, as can be seen in Figure 2 -our method yields lower errors in flow estimation.

Many more experiments -are available in appendix G. We study the impact of key design choices, artificial generalization, and dynamical loss. We show qualitative results on time interpolation, time extrapolation on the Navier dataset. We explore generalization to different grids. We provide more empirical evidence of the soundness of Step 2 in an ablation study (including comparison with attentive neural process Kim et al. ( 2018), an attentionbased structure somehow close to ours), and observe attention maps on several examples. We show that our state estimator goes beyond local interpolation, as conventional interpolation algorithms would do. Finally, we also measure the computational burden of the discussed methods and show that our approach is more efficient.

CONCLUSION

We exploit a double dynamical system formulation for simulating physical phenomena at arbitrary locations in time and space. Our approach comes with theoretical guarantees on existence and accuracy without knowledge of the underlying PDE. Furthermore, our method generalizes to unseen initial conditions and reaches excellent performances outperforming existing methods. Potential applications of our model goes beyond fluid dynamics and can be applied to various PDE-based problem. Yet, our approach relies on several hypotheses such as regular time sampling and observability. Finally, for known and well-studied phenomena, it would be interesting to add physics priors in the system, a nontrivial extension that we leave for future work.

Reproducibility -the detailed model architecture is described in the appendix. For the sake of reproducibility, in the case of acceptance, we will provide the source code for training and evaluating our model, as well as trained model weights. For training, we will provide instructions for setting up the codebase, including installing external dependencies, pre-trained models, and pre-selected hyperparameter configuration. For the evaluation, the code will include evaluation metrics directly comparable to the paper's results.

Ethics statement -While our simulation tool is unlikely to yield unethical results, we are mindful of potential negative applications of improving fluid dynamics simulations, particularly in military contexts. Additionally, we strive to minimizing the carbon footprint associated with our training processes.

A WEBSITE AND INTERACTIVE ONLINE VISUALIZATION

An anonymous website has been created where results can be visualized with an online interactive tool, which allows one to choose time steps interactively with the mouse, and in the case of the Shallow Water dataset, also the orientation of the spherical data:

https://continuous-pde.github.io/

B PROOF OF PROPOSITION 1

The proof proceeds by successive majorations and triangular inequalities. For sake of clarity, and only in this proof we omit the d subscript and write s[n] and z[n] for s d [n] and z d [n], respectively.

We start with ŝ

[n] := ĥ1 • f1 n • ê s[0]
. Thus for any integer n > 0 we have

|s[n] -ŝ[n]| = |h 1 z[n] -ĥ1 ẑ[n] |. ( 14
)
Using Lipschitz property and 4, then

|s[n] -ŝ[n]| ⩽ |h 1 z[n] -ĥ1 z[n] | + | ĥ1 z[n] -ĥ1 ẑ[n] | (15) ⩽ δ h + L h |z[n] -ẑ[n]|.
Noticing that one can rewrite ẑ

[n] as ẑ[n] = f1 n • ê s[0] . Since z[n] = f 1 n z[0]
and using a similar decomposition as for 15), one gets:

z[n] -ẑ[n] ⩽ δ f n-1 k=0 L k f + L n f z[0] -ẑ[0] . (16)
Hence, from equation 15, and using z[0] = e s[0] and ẑ[0] = ê s[0] , we have

s[n] -ŝ[n] ⩽ δ h + L h δ f L n f -1 L f -1 + L n f δ e . ( 17
)
We now move on to the classic auto-regressive case, i.e.

ŝar [n] = ĥ1 • f1 • ê n s[0] . s[n] -ŝar [n] ⩽ h 1 z[n] -ĥ1 z[n] + ĥ1 z[n] -ĥ1 ẑar [n] (18) ⩽ δ h + L h z[n] -ẑar [n] ⩽ δ h + L h δ f + L f e s[n -1] -ê ŝar [n -1] ⩽ δ h + L h δ f + L f δ e + L e s[n -1] -ŝar [n -1] ⩽ δ n-2 i=0 L i + L n-1 s[1] -ŝar [1] , with δ = δ h + L h δ f + L h L f δ e and L = L h L f L e . Moreover, |s[1] -ŝar [1]| = | ĥ1 (z[1]) -f1 ( ẑar [1])| (19) ⩽ δ h + L h |z[1] -ẑar [1]| (20) ⩽ δ h + L h (δ f + L f |z[0] -ẑar [0]|) (21) ⩽ δ (22)
Putting it all together, we get equation 6:

|s[n] -ŝar [n]| ⩽ δ L n -1 L -1 (23)
Finally, equation 17 and equation 23 conclude the proof.

C COMPARISON OF UPPER BOUNDS IN PROPOSITION 1

We start by formulating equation 5 and equation 6 under a comparable form

|s d [n] -ŝd [n]| ⩽ δ + L h L f δ f L n-1 f -1 L f -1 + L h L f δ e (L n-1 f -1) = δ + K 1 (24) |s d [n] -ŝar d [n]| ⩽ δ + δ L n -L L -1 = δ + K 2 (25)
Now we consider two cases depending on the Lipschitz constants of the problem, namely L h , L f , and L e . First, consider the case where the Lipschitz constants are very large (i.e. L h , L f , L e ≫ 1).

In that case, the upper bounds can be approached by

K 1 ≈ L h δ f L n-1 f + L h L n f δ e ( 26
)
K 2 ≈ δ h (L h L f L e ) n-1 + L h δ f L n-1 f L n-1 h L n-1 e + L h L n f δ e L n-1 h L n-1 e (27)
Hence, K 2 ≫ K 1 (we highlighted the difference between both terms in the previous equation. Now consider the case where the Lipschitz constants are very small (i.e. L h , L f , L e ≪ 1). Recall that this case corresponds to a trivial prediction task since any trajectory of System 1 will converge to a unique state. Again, the upper bounds can be approached by

K 1 ≈ 0 (28) K 2 ≈ Lδ (29)
In this trivial case, the upper bound on the prediction error using our method is a combination of the approximation errors from each function. On the other hand, using the classic AR scheme implies a larger error, since the model accumulates approximations at each time step not only from the dynamics but also from the observation function and the encoder.

D PROOF OF PROPOSITION 2

The proof follows the lines of Janny et al. (2022b). The existence of ψ q is granted by the observability assumption. Indeed assumption A2. states that for all q > p, O q is injective in S. Hence, it exists a inverse mapping O * q : O q : S → S such that ∀s

′ ∈ S O * q O q (s ′ ) = s ′ (30) Let z d [0|q]= z d [0] • • • z d [q]
. Hence, one can build ψ q using the dynamics of the system for all x ∈ Ω:

∀s 0 ∈ S, S(s 0 , x, t) = S O * q z d [0|q] , x, t := ψ q z d [0|q], x, t (31)
Now, because of the noise, the disturbed observation ẑd [0|q] = z d [0|q] + δ 0|q may not belong to O q (S), where the inverse mapping O * q is well defined. We solve this by finding the closest "possible" observation. ŝ0 = arg min

s ′ ∈S ẑd [0|q] -O q (s ′ ) (32) ŝ(x, t) = S(ŝ 0 , x, t) := ψ q ẑd [0|q], x, t . (33)
Hence, we have for all

s ′ ∈ S ẑd [0|q] -O q ŝ0 ⩽ ẑd [0|q] -O q s ′ . ( 34
)
In particular, for s ′ = s 0 and since

O q (s 0 ) = z d [0|q], ẑd [0|q] -O q ŝ0 ⩽ ẑd [0|q] -O q s 0 (35)
⩽ δ 0|q . In the other hand, from assumption A2. equation 8: Moreover, since f 2 is Lipschitz

α(p)|ŝ 0 -s 0 | S ⩽ |O q (ŝ 0 ) -O q (s 0 )| (36) ⩽ |O q (ŝ 0 ) -ẑd [0|q]| + | ẑd [0|q] -O q (s 0 )| ⩽ 2 δ 0|q
∂ ∂t |S(s 0 , x, t) -S(ŝ 0 , x, t)| S = |f 2 S(s 0 , x, t) -f 2 S(ŝ 0 , x, t) | S (37)
⩽ L s |S(s 0 , x, t) -S(ŝ 0 , x, t)| S .

and using the Grönwall inequality

|S(s 0 , x, t) -S(ŝ 0 , x, t)| S ⩽ e Lst |s 0 -ŝ0 | S . (38)
Finally, combining equation 36 and equation 38

|S(s 0 , x, t) -S(ŝ 0 , x, t)| S ⩽ 2α(q) -1 δ 0|q e Lst .
which concludes the proof.

E MODEL DESCRIPTION

In this section, we describe the architecture of our implementation in more detail. The model is illustrated in figure 3.

Step 1 -The output predictor derived from System 1 is implemented as a multi-layer graph neural network inspired from [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF]; Sanchez-Gonzalez et al. (2020) but without following the standard "encode-process-decode" setup. Let X = {x 0 , ..., x K } be the set of sub-sampled positions extracted from the known locations X (cf. Artificial generalization from section 3.4). The input of the module is the initial condition at the sampled points and the corresponding positions

x i , sd [0](x i ) i and is encoded into a graph-structured latent space z d [0] = (z d [0] i , e[0] ij ) i,j where z d [0]
i is a latent node embedding for position x i and e[0] ij is an edge embedding for edge pairs (i, j) extracted from a Delaunay triangulation. The encoder ê maps the sparse IC to node and edge embeddings using two MLPs, f edge and f node : Eagle -Eagle is a large-scale fluid dynamics dataset simulating the airflow generated by a drone within a 2D room. We extract sequences of length T = 10 from examples within the dataset, limiting the number of points to 3,000 (vertices were duplicated when the number of nodes fell below this threshold).

z d [0] i = f node sd [0](x i ), x i , e [ 0]ij = f edge x i -x j , |x i -x j | , (39)
The spatially down-sampled versions of these datasets (employed in Table 1 and2) were obtained through masking. We generate a random binary mask, shared across the training, validation, and test sets, to remove a specified number of points based on the desired scenario. Consequently, the observed locations remain consistent across training, validation, and test sets, except Eagle, where the mesh varies between simulations. For Navier and Shallow Water, the High setup retains 25% of the original grid, the Middle setup retains 10%, and the Low setup retains 5%. In the case of Eagle, the High setup preserves 50% of the original mesh, while the Low setup retains only 25%. Temporal down-sampling was also applied by regularly removing a fixed number of frames from the sequences, corresponding to no down-sampling (1/1 setup), half down-sampling (1/2), and quarter down-sampling (1/4). During evaluation, the models are tasked with predicting the solution to every location and time instant present in the original simulation.

G MORE RESULTS

Time continuity -is illustrated in Figure 6 on the Navier dataset. Our model and the baselines are trained in a very challenging setup, where only part of the information is available. During training, not only does the spatial mesh only contains 25% of the complete simulation grid, but also the time-step is increased to four time its initial value. In this situation, the model needs to represent low-resolution data while being trained on sparse data.

Generalization to unseen future timesteps -Beyond time continuity, our model offers some generalization to future timesteps. Table 3 shows extrapolation results for high/mid/low subsampling of the spatial data on the Navier dataset which outperforms the predictions of competing baselines.

Generalization to unseen grid -In our spatial and temporal interpolation experiments (tables 1 and 2 of the main paper), we assumed that the observed mesh remains identical during training and testing. Nevertheless, the ability to adapt to diverse meshes is an important aspect of the task. To evaluate this capability, we trained our model in the spatial extrapolation setup on the Navier dataset. We compute the error when exposed to different meshes, potentially with a different sampling rate, and report the results in table 4. Our model demonstrates good generalization skills when confronted with new and unseen grids. We observe that the error on new grids is close to the error reported in table 1 in the Ext-X case, we show additionally that the model can generalize even if the observed grid is different. Notably, the model performs well when trained with a medium sampling rate. Despite some performance degradation when the evaluation setup is significantly different compared to training, our model effectively maintains its interpolation quality between out-of-domain error (Ext-X ) and in-domain error, testifying to the robustness of our dynamic interpolation module.

Ablations -we study the impact of key design choices in Figure 7a. First, we show the effect of the subsampling strategy to favor learning of spatial generalization, c.f. Section 3.4, where we subsample the input to the auto-regressive backbone by keeping 75% of the mesh. We ablate this feature by training the model on 100%, 50%, and 25% of the input points. When the model is trained on 100% of the mesh, it fails to generalize to unseen locations, as the model is always queried on points lying in the input mesh. However, reducing the number of input points significantly further from the operating point decreases the performance of the backbone, as it does not dispose of enough points to learn meaningful information for prediction. We also replace the final GRU with simpler aggregation techniques, such as a mean and a maximum pooling, which drastically degrades the results. Finally, we ablate the dynamics part of the training loss (Eq. 13). As expected, this deteriorates the results significantly.

More ablation on the interpolator -We conducted an ablation study to show that limiting attention is detrimental. To do so, we designed four variants of our interpolation module:

• Single attention (w/o GRU) -performs the attention between the query and the embeddings in a single shot, rather than time-step per time-step. This variant neglects the insights from control theory presented in section 3.1 (Step 2). The single softmax function limits the attention to a handful of points, whereas our method encourages the model to attend to at least one point per time step and reason on a larger timescale, considering past and future predictions, which is beneficial for interpolation tasks, as supported by proposition 2.

• Spatial (w/ GRU) & Temporal (w/o GRU) neighborhood -limit the attention to the nearest temporal or spatial points, which significantly degrades the metrics. To handle setups with sparse and subsampled trajectories, the interpolation module greatly benefits from not only distant points but also from the temporal flow of the simulation. • Temporal attention (w/o GRU) replaces the GRU in our model with a 2-head attention layer. This variant of our model does not improve the performance compared to a GRU. We argue that GRU is more suited for accumulating observations in time, as its structure matches classic observer designs in control theory.

In-T , In-X In-T , Ext-X Ext-T , In-X Ext-T , Ext-X (
• Attentive Neural Process Kim et al. ( 2018) is a interpolation module close to ours resembling the Single attention ablation, with an additional global latent c to account for uncertainties. The model involves a prior function q(c, s) trained to minimize the Kullback-Leibler divergence between q z, s X , T (computed using the physical state at observed points) and q c, s Ω \ X , 0, T (computed using the physical state at query points).

Results are shown in table 5. All ablations exhibit worse performance than ours. Note that the ANP ablation involves performing the interpolation in the physical space to compute the Kullback-Leibler divergence during training. Thus, the interpolation module cannot use the latent space from the autoregressive module, which may explain the drop in performance. Adaptating ANP to directly leverage the latent states is probably possible, but not straightforward and requires several key changes in the architecture.

Efficiency -the design choices we made led to a computationally efficient model, compared to prior work. For all three baselines, the required number of computed time steps for the auto-regressive rollout depends on (1) the number of predicted time steps, and (2) the time values themselves, as for later values of t, more iterations need to be computed. In contrast, our method forecasts using attention from a set of "anchor states", which is controlled through the hyper-parameter ∆. The length of the auto-regressive rollout is therefore constant and does not depend on the number of predicted time steps. Furthermore, while DINo scales very well to predict additional locations, it requires a costly optimization step to compute α 0 . MGN does benefit from the efficient cubic interpolation algorithm, which is a side effect of the fact that it has been adapted to this task, but not

t0 t0 + ∆ t0 + 2∆ t0 + 3∆ t0 + 4∆ t0 + 5∆ t0 + 6∆ 0.002 0.004 ∂ ŝ(x,τ ) ∂zd[n] 2 (a)
Frontier tracking: when queried on a streamline between areas of opposite vorticity, the interpolation module attends not only to the spatial neighborhood but also to the temporal flow near the frontier. designed for it. We experimentally confirm these claims in Figure 7, where we provide the evolution of runtime as a function of query locations, and of query time steps, respectively. In both cases, our model compares very favorably to competing methods.

Attention maps -To further support our claims, we analyzed the behavior of the interpolation module in more depth and showed the top-100 most important nodes from the embedding points z d [n](x i ) used to interpolate at different queries. The figure is shown in Figure 8. We observed very complex behaviors that dynamically adapt to the global situation around the queried points.

Our interpolation module appears to give more importance to the flow rather than merely averaging the neighboring nodes, thus relying on "why" the queried point is in a specific state. Such behavior would be extremely difficult to implement in a handcrafted algorithm.

Parameter Sensitivity Analysis -We investigate the influence of two principal hyper-parameters, namely the step size ∆ and the number of residual GNN layers L, on the performances of our model. We present the results of our experiments in figure 9 on the Navier dataset, which has been spatially down-sampled at 10% during training and has a temporal resolution reduced by two.

The choice of the step size between iterations of the auto-regressive backbone directly affects both training and inference time. For a trajectory of T frames, the number of anchor states z d [n] is determined by ⌊T /∆⌋. Increasing the step size ∆ of the learned dynamics leads to a higher number of embeddings over which the models need to reason. A parallel can be drawn between this phenomenon and the influence of the discretization size on the accuracy of numerical methods for solving PDEs. Furthermore, the selection of ∆ also impacts the generalization capabilities of the model in Ext-T . When ∆ > ∆ * , the model is queried during training on intermediate instants not directly associated with any of the anchor states z d [n]. This is visible in Figure 9 where, for instance, with ∆ = ∆ * , the In-X /In-T error is the lowest, but other metrics increases compared to ∆ = 2∆ * .

The number of layers L in the auto-regressive backbone significantly influences the overall performance of the model, both within the domain and on the exteriors. Increasing the number of layers generally leads to improved performance. However, it appears that beyond L = 8, the error starts to Figure 9: Impact of hyper-parameters on model performance -We evaluate the impact of three critical hyper-parameters on our architecture, namely, (a) the step size ∆, the depth L of the autoregressive backbone and (b) the weighting of the dynamics cost in equation 13. To assess the performance, we employ the 10% Navier dataset with 1/2 frames and compute metrics for both in-domain and out-domain. The results reveal that increasing the depth of the GNN layers enhances the model's performance, while lower values of ∆ lead to better metrics. However, we observed a degradation in the ability of the model to generalize to unseen time instants for the special case ∆=∆ * . Moreover, we found that equally weighting of both terms L continuous and L dynamics leads to best results. increase, indicating a saturation point in terms of performance gain. The relationship between the number of layers and model performance is visually depicted in Figure 9. Throughout this paper, we maintain this hyper-parameter constant for the sake of simplicity, as our primary focus is the spatial and temporal generalization of the solution.

Failure cases -we expose failure cases on the Eagle dataset (in the Low spatial down-sampling scenario) in figure 10. In some particularly challenging instances of this turbulent dataset, we noticed drops in accuracy located in fast-evolving regions of the simulation, and in particular near the flow source. We hypothesize that the origin of the failure is related to the comparatively smaller processor unit used in our auto-regressive backbone compared to the baseline introduced in Janny et al. ( 2023), hence producing less accurate anchor states when the horizon increases.

Figure 10: Failure cases on Eagle-We observed failure cases on highly challenging instances of the EAGLE dataset as the prediction horizon is increasing. We show the per point error in three different instances and observed that the error increases with the time horizon, especially close to turbulent areas, such as below the UAV.

…Figure 1 :

 1 Figure1: Model overview -We achieve space and time continuous simulations of physics systems by formulating the task as a double observation problem. System 1 is a discrete dynamical model used to compute a sequence of latent anchor states z d auto-regressively, and System 2 is used to design a state estimator ψ q retrieving the dense physical state at arbitrary locations (x, t).

Figure 2 :

 2 Figure 2: Results on Eagle -Per point error of the flow prediction on an Eagle example in the Low spatial down-sampling scenario. Our model exhibits lower errors as also shown in Tables1 and 2.

eFigure 3 :

 3 Figure3: Model overview -The model leverages a dynamical system (System 1) to perform autoregressive predictions of the dynamics in a mesh-structured latent space from sparse initial conditions. It is combined with a data-driven state estimator derived from another continuous-time dynamical system (System 2), implemented with multi-head cross-attention. The attention mechanism queries the intermediate anchor states from the auto-regressive predictor and uses Fourier positional encoding to encode the query points (x, τ ). An additional GRU refines the dynamics after interpolation.

Figure 4 :

 4 Figure4: MaGNet -suffers from drastic shifts in distribution between training and evaluation. The model is trained on points from X , which corresponds to a small portion of the domain. We used our subsampling trick to artificially generate queries. During evaluation, we require the prediction at every available point in the complete simulation, hence, MaGNet must interpolate the initial condition to a large number of query points, filling the input of the auto-regressive model with noisy estimates of the IC.

Figure 5 :

 5 Figure 5: Qualitative results on Shallow-Water -Simulation obtained with our model and the baseline in the challenging 5% setup on the Shallow Water dataset (without temporal sub-sampling). Each model is initialized with a small set of sparse observations and needs to extrapolate the solution at many unseen positions. Our model outperforms the baselines, which struggle to compute the solution outside the training domain.

Figure 6 :

 6 Figure 6: Time continuity on the Navier dataset -during training, models are only exposed to a sparse observation of the trajectories, represented spatially by the dots in the upper left figure and temporally by the semi-transparent frames. Our model maintains the temporal consistency of the solution and outperforms the baselines.

Figure 7 :

 7 Figure 7: Ablations and runtime -(a) Ablations on Navier (Yin et al., 2022; Stokes, 2009) with 10% of data and half temporal resolution, from left to right: exploring subsampling strategies, replacing GRU par mean/max pooling, removing physics grounding. (b) Runtime analysis as a function of query locations and time steps, respectively. The graph shows the average runtime (over 100 runs), and shadows indicate lower and upper bounds over the runs.

  Blob tracking: in homogeneous areas, the model tracks the origin of the perturbation, and focuses on its displacement. Our dynamic interpolation exploits the evolution of the state rather than merely averaging neighboring nodes. Periodic boundaries: our model effectively leverages the periodic condition of the Navier dataset, especially when queried on points originating from perturbations on the other side of the simulation. Again, the interpolation depends on which points explain the output, rather than the neighborhood.

Figure 8 :

 8 Figure 8: Norm of output derivative -wrt. each z d [n](x i ) (Navier, high spatial subsampling setup). We display top-100 nodes (•) with the highest norm, i.e. the most important nodes for interpolation at query point (■). Using gradients rather than attention allows us to visualize the action of the GRU. We observe context-adaptive behaviors, leveraging temporal flow information over local neighbors, challenging to implement in handcrafted algorithms.

Table 1 :

 1 is a standard multi-layered GNN Space Continuity -we evaluate the spatial interpolation power of our method vs. the baselines and standard interpolation techniques. We vary the number of available measurement points in the data for training from High (25% of simulation grid), Middle (10%), and Low (5%) amount of points and show that our model outperforms the baselines. Evaluation is conducted over 20 frames in the future (10 for Eagle) and we report the MSE to the ground truth solution (×10

				Navier		Shallow Water	Eagle
			High	Mid	Low	High	Mid	Low	High		Low
	DINo	In-X	1.557	1.130	1.878 0.1750 0.1814 0.2733 287.3		302.7
	(Yin et al., 2022)	Ext-X	1.600	1.253	5.493	4.638	13.40	21.55	381.7		489.6
	Interp. MGN	In-X	1.913 0.9969 0.6012 0.3663 0.2835 0.7309 64.44		83.58
	(Pfaff et al., 2020)	Ext-X	2.694	4.784	14.80	1.744	4.221	8.187	173.4		241.5
	Time Oracle (n.c)	In-X Ext-X	0.851	n/a 4.204	15.63	1.617	n/a 4.327	8.522	147.0	n/a	221.2
	MAgNet	In-X	18.17	6.047	8.679 0.3196 0.3358 0.4292 99.79		124.5
	(Boussif et al., 2022) Ext-X	35.73	26.24	57.21	10.21	23.20	30.55	194.3		260.7
	Ours	In-X Ext-X	0.1989 0.2136 0.2446 0.2940 0.3139 0.2700 70.02 0.2029 0.2463 0.5601 0.4493 1.051 2.800 90.88		78.83 117.2

Table 3 :

 3 Time Extrapolation -We assessed the performances of our model vs. the baselines in a time-extrapolation scenario by forecasting the solution on a horizon two times longer than the training one (i.e. 40 frames). Our model remains more performant.

						Training
					Navier		Shallow
				High	Mid	Low	High	Mid	Low
		In-X High Ext-X	0.2492 0.7929 4.5165 0.5224 1.5431 4.3447 0.2477 0.7782 4.4038 0.5256 1.5822 4.4963
	Evaluation	Mid	In-X Ext-X	0.4370 0.3230 0.9759 0.8528 1.2908 3.6766 0.4410 0.3401 0.9496 0.8617 1.2589 3.6043
		Low	In-X Ext-X	2.2000 0.4039 0.6732 2.4395 1.5634 3.4793 2.2037 0.4216 0.7892 2.3914 1.5313 3.2334

Table 4 :

 4 Generalization to unseen grid -We investigate generalization to previously unseen grids by training our model on the Navier dataset in the space extrapolation setup. We report the error (MSE (×10 -3 )) inside and outside the spatial domain X measured with different sampling rates unseen during training. The diagonal shows results on grids with identical sampling rates wrt. training, but sampled differently. Our model shows great generalization properties.

Table 5 :

 5 Ablation on interpolation -We performed four ablations on the interpolation module (MSE (×10 -3 )). Single attention combines all z d [n] into a single key vector, employing attention only once (w/o GRU). Temporal attention replaces the GRU with a 2-head attention, Spatial neigh. restricts attention to the five spatially nearest points from the query, and Temporal neigh. computes attention only with the nearest time z d [n] to the queried time τ (w/o GRU). These results indicate that considering long-range spatial and temporal interactions is beneficial for the interpolation task.

						Ablation Study
	Mean Squared Error	10 -3
		10 -4	100%	50%	25%	Mean	Max	w/o grounding.	Ours
				Spatial Sub.		Agregation

In what follows, while being a misnomer, space and time continuity of the solution designate the continuity of the spatial and temporal domain of definition of the solution, and not the continuity of the solution itself.

Code will be made public. Project page: https://continuous-pde.github.io/

Since the simulation is conducted up to T , and considering the time step ∆, in practice q ⩽ ⌊ T ∆ ⌋

ACKNOWLEDGEMENTS

We recognize support through French grants "Delicio" (ANR-19-CE23-0006) of call CE23 "Intelligence Artificielle" and "Remember" (ANR-20-CHIA0018), of call "Chaires IA hors centres". This work was performed using HPC resources from GENCI-IDRIS (Grant 2023-AD010614014).

f nodes and f edges are two ReLU-activated MLPs, each consisting of 2 layers with 128 neurons. The initial node and edge features z d [0] i and e[0] ij are represented as 128-dimensional vectors.

The dynamics f1 is modeled as a multi-layered graph neural network inspired from [START_REF] Pfaff | Learning mesh-based simulation with graph networks[END_REF]; [START_REF] Sanchez-Gonzalez | Learning to simulate complex physics with graph networks[END_REF], we therefore add a layer superscript ℓ to the notation:

The GNNs employ two MLPs g ℓ node and g ℓ edges with same dimensions as f edges and f nodes . We compute the sequence of anchor states

in the latent space by applying f1 auto-regressively.

The observation function ĥ1 extracts the sparse observations sd [n] from the latent state z d [n] and consists of a two-layered MLP with 128 neurons, with Swish activation functions [START_REF] Ramachandran | Searching for activation functions[END_REF] applied on the node features, i.e. sd [n]

Step 2 -The spatial and temporal domains Ω × 0, T are normalized, since it tends to improve generalization on unseen locations. The state estimator ψ q takes as input the sequence of latent graph representation

and a spatiotemporal query sampled in Ω × 0, T . This query is embedded in a Fourier space using the function ζ ω which depends on a frequency parameter ω ∈ R dim Ω+1 (initialized uniformly in [0, 1]). By concatenating harmonics of this frequency up to some rank, we obtain a resulting embedding of 128 dimensions (if ζ ω (x, t) exceeds the number of dimensions, cropping is performed to match the target shape).

(41) The continuous variables z n∆ (x, t) conditioned by the anchor states are computed with a multi-head attention [START_REF] Vaswani | Attention is all you need[END_REF]

where f mha is defined as

Here, A(•, •, •) refers to the multi-head attention mechanism described in [START_REF] Vaswani | Attention is all you need[END_REF] with four attention heads, and B(•) represents a single-layer multi-layer perceptron activated by the rectified linear unit (ReLU) function. We do not use layer normalization.

The Gated Recurrent Unit [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF] aggregates the sequence of conditioned variables (of length q) as follows:

where u[n] is the hidden memory of a GRU, initialized at zero. r GRU denotes the update equations of a GRU -we omit gating functions from the notation -and D is a decoder MLP that maps the final GRU hidden state to the desired output, that is, the value of the solution at the desired spatiotemporal coordinate (x, t), We used a two-layered gated recurrent unit with a hidden vector of size 128, and a two-layered MLP with 128 neurons activated by the Swish function for D.

Training loop -To create artificial generalization scenarios during training, we employ spatial subsampling. Specifically, during each gradient iteration, we randomly and uniformly mask 25% of X and feed the remaining 75% to the output predictor (System 1). To reduce training time further and improve generalization on unseen locations, we use bootstrapping by randomly sampling a smaller set of points for querying the model (i.e. as inputs to ψ q ). To do so, we maintain a probability weight vector W of dimension |X ×T |, initialized to one. At each gradient descent step, we randomly select N =1, 024 points from X × T weighted by W . We update the weight matrix by setting the values at the sampled locations to zero and then adding the loss function value to the entire vector. This procedure serves two purposes: (a) it keeps track of poorly performing points (with higher loss) and (b) it increases the sampling probability for points that have been infrequently selected in previous steps.

The choice of ∆ in the dynamics loss equation 13 allows us to reduce the complexity of the model.

In Table 1, we present results obtained with ∆ = 3∆ * indicating that the output predictor (System 1) predicts the latent state representation three time steps later. Consequently, the number of autoregressive steps during training decreases from T /∆ * (e.g., for MeshGraphNet and MAgNet) to T /∆. In Table 2, we used ∆ = 2∆ * . For a more comprehensive discussion on the effect of ∆ on performances, please refer to Appendix G.

Training parameters -To be consistent, we trained our model with the same training setup over all different experiments (i.e. same loss function, and same hyper-parameters). However, for the baseline experiments, we did adapt hyper-parameters and used the ones provided by the original authors when possible (see further below). We used the AdamW optimizer with an initial learning rate of 10 -3 . Models were trained for 4,500 epochs, with a scheduled learning rate decay multiplied by 0.5 after 2,500; 3,000; 3,500; and 4,000 epochs. Applying gradient clipping to a value of 1 effectively prevented catastrophic spiking during training. The batch size was set to 16.

F BASELINES AND DATASETS DETAILS F.1 BASELINES

The baselines are trained with the AdamW optimizer with a learning rate set at 10 -3 for 10,000 epochs on each dataset. We keep the best-performing parameters on the validation set for evaluation on the test set.

DINo -we used the official implementation and kept the hyper-parameters suggested by the authors for Navier and Shallow Water. For Eagle, we used the same hyper-parameters as for Shallow Water.

The training procedure was left unchanged.

MeshGraphNet -we used our own implementation of the model in PyTorch, with 8 layers of GNNs for Navier and Shallow Water, and up to 15 for Eagle. Other hyper-parameters were kept unchanged. We warmed up the model with single-step auto-regressive training with noise injection (Gaussian noise with a standard deviation of 10 -4 ), as suggested in the original paper, and then fine-tuned the parameters by training on the complete available horizon. Both steps try to minimize the mean squared error between the prediction and the ground truth. Edges are computed using Delaunay triangulation. During evaluation, we perform cubic interpolation between time steps (linear interpolation gives better results on Eagle) first, then 2D cubic interpolation on space to retrieve the complete mesh.

MAgNet -We used our own implementation of the MAgNet[GNN] variant of the model, and followed the same training procedure as for MeshGrapNet. The parent mesh and the query points are extracted from the input data using the same spatial sub-sampling technique as in ours, and the edges are also computed with Delaunay triangulation. During evaluation, we split the query points into chunks of 10 nodes, and compute their representation with all the available measurement points. This reduces the number of interpolated vertices in the input mesh and improves performances at the cost of higher computation time (see figure 4). However, to be fair, this increase in computational complexity introduced by ourselves was not taken into account when we discussed computational complexity in appendix G.

F.2 DATASET DETAILS

Navier & Shallow Water -Both datasets are derived from the ones used in [START_REF] Yuan Yin | Continuous pde dynamics forecasting with implicit neural representations[END_REF]. We adopted the same experimental setup but generated distinct training, validation, and testing sets. For details on the GT simulation pipeline, please see [START_REF] Yuan Yin | Continuous pde dynamics forecasting with implicit neural representations[END_REF]. The Navier dataset comprises 256 training simulations of 40 frames each, with additional two times 64 simulations allocated for validation and testing. Simulations are conducted on a uniform grid of 64 by 64 pixels (i.e. Ω),