Quentin Goestchel

Gwenaël Guillaume
email: gwenael.guillaume@cerema.fr

David Ecotière

Benoit Gauvreau

Analysis of the numerical properties of the transmission line matrix model for outdoor sound propagation

Keywords: Transmission Line Matrix Method, Numerical Methods, Outdoor Sound Propagation, Dispersion

Outdoor sound propagation modeling can be performed through various numerical models. Nowadays, time-domain methods are widely used and enable accurate simulations for most of the sound physical phenomena involved in environmental acoustics. However, the transmission line matrix (TLM) method remains relatively poorly documented. In this paper, a thorough review of the method and a stability analysis are provided. The review of the model leads to a more robust understanding of its link with the wave equation using Taylor expansions. Two different cases are considered: a homogeneous non-dissipative medium and an inhomogeneous dissipative one. The stability analysis shows the quantification of the model inherent dispersion error and highlights similarities with a finite difference scheme. To conclude, a numerical experiment aims to characterize precisely the limitations of the method applied to outdoor sound propagation.

Introduction

Anthropogenic noise pollution is known to have a negative impact on public health and on biodiversity [START_REF]Burden of Disease from Environmental Noise: Quantification of Healthy Life Years Lost in Europe[END_REF][START_REF] Warren | Urban bioacoustics: It's not just noise[END_REF][START_REF] Sordello | Evidence of the impact of noise pollution on biodiversity: A systematic map[END_REF]. Disturbing sounds generated by human activities are not confined to urban environments and reach further into natural areas. Predictions of long-range outdoor sound propagation are useful for assessing the impact of noise pollution on living beings and their environment. For this purpose, numerical acoustic models are necessary and a large number of them are available in the literature [START_REF] Hornikx | Ten questions concerning computational urban acoustics[END_REF]5]. From frequency domain methods to time domain methods, the study and improvement of these models have recently developed with the available computing power.

Unlike other time-domain numerical methods such as finite differences or the Fourier pseudospectral time-domain method, the transmission line matrix method (TLM) has not been documented extensively [5][START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF][START_REF] Hornikx | openPSTD: The open source pseudospectral time-domain method for acoustic propagation[END_REF]. Originally introduced to model high frequency electromagnetic fields, the TLM model and its formalism were introduced in 1971 [START_REF] Johns | Numerical solution of 2-dimensional scattering problems using a transmission-line matrix[END_REF][START_REF] Johns | On the Relationship Between TLM and Finite-Difference Methods for Maxwell's Equations (Short Paper)[END_REF]. In his book, de Cogan showed the flexibility of this method by applying it to various domains and for electromagnetic waves, the TLM model seems to tackle with most of the difficulties usually encountered by numerical methods such as perfectly matched layers (PML), moving boundaries, or non-cartesian meshes [START_REF] De Cogan | Transmission Line Matrix (TLM) in Computational Mechanics[END_REF][START_REF] Guillaume | Time-domain impedance formulation for transmission line matrix modelling of outdoor sound propagation[END_REF].

In the field of acoustics, the TLM model is successfully used in room acoustics and the essence of the method has been studied from a theoretical point of view [START_REF] Kagawa | Discrete Huygen's model approach to sound wave propagation[END_REF][START_REF] Kagawa | Discrete Huygens' model approach to sound wave propagation-reverberation in a room, sound source identification and tomography in time reversal[END_REF]. For outdoor scenarios, application to inhomogeneous media, urban areas or even forest-like environments have been developed and compared to experimental results [START_REF] Kagawa | Discrete Huygens' modelling simulation of sound wave propagation in velocity varying environments[END_REF][START_REF] Guillaume | Use of the transmission line matrix method for the sound propagation modeling in urban area[END_REF][START_REF] Chobeau | A Transmission Line Matrix model for sound propagation in arrays of cylinders normal to an impedance plane[END_REF]. The numerical scheme has even been optimized by making it capable of running on parallel graphic processing units [START_REF] Guillaume | Optimized transmission line matrix model implementation for graphics processing units computing in built-up environment[END_REF]. Finally, complex features such as impedance boundary conditions or meteorological effects have been introduced into the TLM model in order to model more realistic environments [START_REF] Guillaume | Time-domain impedance formulation for transmission line matrix modelling of outdoor sound propagation[END_REF][START_REF] Aumond | Application of the Transmission Line Matrix method for outdoor sound propagation modelling -Part 2: Experimental validation using meteorological data derived from the meso-scale model Meso-NH[END_REF]. However, no stability analysis of the method has been undertaken yet. Therefore, this paper aims to quantitatively define the limitations of the TLM method through an extensive theoretical analysis. The scope of the study focuses on long-range outdoor noise predictions, as they are complex situations to model and therefore the most likely to expose limitations of numerical models. For this purpose, only the central scheme is considered, leaving aside boundary conditions, PMLs and other artifacts inducing numerical errors.

An updated theoretical d-dimension formalism of the TLM is first proposed in the next section. In the third and fourth section, the model is detailed for a homogeneous non-dissipative and an inhomogeneous dissipative case, respectively. Within each part, the link with the wave equation is demonstrated and a stability analysis is performed, highlighting the inherent dispersion error of the model. In a last section, a numerical experiment is presented and an analysis of the results is given and interpreted in order to summarize precisely the limitations of the TLM model use for outdoor sound propagation.

Theoretical basis

Overall formalism

Unlike most numerical methods, the TLM is not initially derived from the conservative first order equations. It is based on the Huygens' principle for wave propagation and on an electro-acoustical analogy [START_REF] Kagawa | Discrete Huygen's model approach to sound wave propagation[END_REF]. The property of a wavefront to be recursively discretized in the form of punctual secondary sources is exploited and applied to a Cartesian mesh. The sound propagation is then modeled as pressure pulses, propagating along transmission lines. This decomposition of the propagation mechanism is the TLM basis. It allows the analogy between the progression of a sound wave and the diffusion of pulses between the nodes of a mesh.

In the following, a d-dimension spatial generalization inspired from previous works is described [START_REF] Guillaume | Application of the transmission line matrix method for outdoor sound propagation modelling -Part 1: Model presentation and evaluation[END_REF]. The spatial domain is equally discretized in every directions with steps ∆ such as x d = j d ∆ , j d ∈ Z and the time is decomposed into steps such as t n = n∆t, n ∈ N. The spatial locations of nodes on the grid are then represented by the vector of indices: r = (j 1 , . . . , j d ). The core of the model is that each node in the volume under consideration receives and emits incident and scattered pulses instantaneously, at each time step n. These pressure pulses are travelling along links between the nodes, named transmission lines m. The formalism for this representation is given as follow: n I m r is the incident pressure pulse to the node r along the line m at the time step n and n S m r is the reflected instantaneous scattered pulse.

Scattering matrix

Expressing this mechanism in d-dimensions leads to the main asset of the TLM method: the scattered pulses are calculated only with the matrix relation:

n S r = n D r • n I r , (1)
with n S r and n I r the scattered and incident pulses vectors and n D r the diffusion matrix

n D r = n          R 1 r T 2 r • • • • • • T M r T 1 r R 2 r T m r • • • . . . . . . T 2 r . . . T m r . . . . . . . . . T m r R m r T M r T 1 r • • • • • • T m r R M r          r , (2)
with m ∈ M , the number of lines around a node. The M × M matrix is populated with the pressure reflection and transmission coefficients, representing the behavior of pressures pulses encountering impedance discontinuities at a node, leading to reflections in the incident transmission lines and "scattering" in the others, whose impedance coefficients are formulated as Z m Lr . Z m Tr is the impedance of the equivalent medium after the node discontinuity, as illustrated by Fig. 1(a). Fig. 1 represents these mechanisms in the 2D case of an incident pulse to a node from a single transmission line. The general expressions of the coefficients are:

R m r = Z m Tr -Z m Lr Z m Tr + Z m Lr , (3a)
T m r =1 + R m r . (3b)

Connection laws

To complete the method, connection laws, describing the time relation between the scattered pulses at t n and incident pulses at t n+1 are needed. A d-dimension generalization of these laws can be written as:

70 for k ∈ {1, . . . , d} : n+1 I 2k-1 (j1,...,j d ) = n S 2k (j1-δ k1 ,...,j d -δ kd ) n+1 I 2k (j1,...,j d ) = n S 2k-1 (j1+δ k1 ,...,j d +δ kd ) , (4)
δ being used here as the Kronecker delta. A 2D example of these laws is illustrated by Fig. 2.

Top node pulse

Figure 2: Representation of the 2D connection laws at the local node r = (j 1 , j 2 ). Scattered pulses "S" (orange) at time step n become incident pulses "I" (blue) at step n + 1.

Eqs. (1) to (4) are essentials to describing and implementing the TLM model. The link between them and the wave equation will be demonstrated in order to legitimize the resolution method. To do so, a general version of Eq. (1) for each 75 line is needed; it can be expressed as:

n S m r = n R m r n I m r + M m=1 n T m r n I m r -T m r n I m r , m ∈ {1, . . . , M }. (5)
Introducing the discrete pressure at the nodes n P r = M m=1 T m r n I m r leads to:

n S m r = n P r -n I m r , ∀m ∈ {1, . . . , M }. (6)
The following sections detail the numerical scheme in the cases of a homogeneous and an inhomogeneous medium and show that Eq. ( 6) and the discrete pressure expression are valid for both.

Scheme for homogeneous non-dissipative media

In this case, the number of transmission lines m ∈ {1, . . . , M = 2×d}, as two nodes are linked to the local node in each spatial direction. Thanks to an electroacoustical analogy [START_REF] Scott | An improved transmission line matrix model for the 2D ideal wedge benchmark problem[END_REF] illustrated on Fig. 3, the reflection and transmission coefficients of the pressure at each node can be expressed by inserting

Z m Lr = Z and Z m Tr = Z 2d -1 in Eq. (3)
:

R m r = 1 -d d , (7a)
T m r = 1 d . ( 7b
)
Figure 3: Electrical analogy scheme along one transmission line in the case of a homogeneous non-dissipative medium.

Link with the wave equation

Using the Millman's theorem [START_REF] De Cogan | Transmission Line Matrix (TLM) in Computational Mechanics[END_REF], the discrete pressure can be expressed as:

n P r = 1 d 2d m=1 n I m r . (8)
Writing the expression (8) at the the time step n + 1 and applying twice the Eqs. ( 4) and ( 6) leads to [START_REF] Hofmann | Simulation of outdoor sound propagation with a transmission line matrix method[END_REF]:

n+1 P r + n-1 P r = 1 d d m=1 n P (j1+δm1,...,j d +δ md ) + n P (j1-δm1,...,j d -δ md ) . (9)
In order to know the approximation order of this scheme, and to retrieve the wave equation, Taylor expansions can be used [START_REF] Kagawa | Discrete Huygen's model approach to sound wave propagation[END_REF]. It is important to point out that n P (j1,...,j d ) is an approximated value of the exact pressure p (x 1 , . . . , x d , t n ) taken at the point (x 1 , . . . , x d ) at time t n . Eq. ( 9) becomes:

∂ 2 p ∂t 2 -c 2 TLM ∇ 2 p = O ∆t 2 + O ∆ 4 ∆t 2 , (10)
with

c TLM = ∆ √ d∆t .
It is now possible to recognize the wave equation and to observe that the TLM model is a second-order approximation method in time, and space. Another interesting observation is that the TLM solves the wave equation only if the condition c 2 TLM = c 2 0 is fulfilled. Developing this condition gives:

c 0 ∆t ∆ = 1 √ d , (11)
which corresponds with the CFL criterion of the so-called finite difference Leap-Frog scheme [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF].

Stability analysis

To evaluate the model stability, the numerical dispersion relation of the method can be written by inserting the formulation of a discretized plane wave

n P r = A exp i d l=1 k l j l ∆ -ωn∆t
in Eq. ( 9). Namely:

cos(ω∆t) = 1 d d m=1 cos(k m ∆ ), ∀∆t, ∀∆ . (12)
This equation implies that the TLM method is unconditionally stable in the homogeneous non-dissipative case. However, this relation also shows that the model presents numerical dispersion in the main directions of the grid (for instance in the 2D case θ = β π 2 , β ∈ Z). To illustrate this phenomenon, the dispersion error is represented as a function of the angle θ between the plane wave vector and the horizontal direction of the mesh. Fig. 4 displays three different cases of spatial grid discretizations that are defined as the numbers of nodes per wavelength along the Cartesian directions. These plots show that the relative dispersion error is below 5% when the number of points per wavelength N is greater than 5.

Scheme for inhomogeneous dissipative media

In order to adapt the model to realistic outdoor scenarios, the inhomogeneous sound velocity and the atmospheric attenuation are modeled by adding two specific transmission lines [START_REF] Kagawa | Discrete Huygens' modelling simulation of sound wave propagation in velocity varying environments[END_REF][START_REF] Guillaume | Application of the transmission line matrix method for outdoor sound propagation modelling -Part 1: Model presentation and evaluation[END_REF]. Their impedance can vary in both space and time thanks to additional parameters η and ζ defined later in the text (see subsection 4.2 to subsection 4.4). A 2D example of the environment around one node is exposed on Fig. 5. The main differences with the homogeneous and non-dissipative medium case (section 3) lie in the reflection and transmission coefficients in the diffusion matrix n D r and the introduction of an additional connection law that are successively detailled in the three following paragraphs.

Scheme along a regular line

The lines whose superscripts m ∈ {1, . . . , 2d} behave identically as those presented in section 3. Their impedance depends only on the medium to model.

However, as exposed on Fig. 6, the appearance of two additional lines changes the equivalent impedance encountered by the pulses at the corresponding nodes and leads to the two following expressions of the diffusion matrix coefficients:

n R m r = - n η r + n ζ r + 2 (d -1) n η r + n ζ r + 2d , ( 13a
) n T m r = 2 n η r + n ζ r + 2d . ( 13b
)

Scheme along a modified impedance line

In order to model local heterogeneity in the sound propagation velocity such as induced by wind or temperature vertical gradients, a transmission line with a reflective termination of length ∆ 2 and with superscript m = 2d + 1, is intro-130 duced. The impedance of this branch is set as

Z 2d+1 r = Z η
, with η a parameter related to the local effective velocity [START_REF] Guillaume | Application of the transmission line matrix method for outdoor sound propagation modelling -Part 1: Model presentation and evaluation[END_REF]. As illustrated on Fig. 5, this line is not linked to any other node and is introduced to simulate a volume variation at the node leading to the following modified expressions of the diffusion matrix coefficients [START_REF] Kagawa | Discrete Huygens' modelling simulation of sound wave propagation in velocity varying environments[END_REF]:

135 n R 2d+1 r = - n η r -n ζ r -2d n η r + n ζ r + 2d , ( 14a
) n T 2d+1 r = 2 × n η r n η r + n ζ r + 2d . ( 14b
)
This extra branch at the node implies an additional connection law to the homogeneous and non-dissipative medium case (Eq. ( 4)) that differs slightly from the usual ones since the scattered pulses are reflected back to the node of origin of the pulse, i.e:

n+1 I 2d+1 r = n S 2d+1 r . (15)

Scheme along an anechoic transmission line

One method to include atmospheric absorption in the TLM model without increasing significantly computational costs consists in inserting an anechoic branch of superscript m = 2d + 2 and impedance Z 2d+2 dependent. This is the first limitation of this method because it assumes a constant absorption coefficient for the frequency range of the simulation.

There is no reflection nor transmission on this line because its purpose is to simulate an amplitude loss of the wave. It is then straightforward to express the reflection and the transmission coefficients as:

n R 2d+2 r = n T 2d+2 r = 0. ( 16
)

Link with the wave equation

Thanks to the previous developments, the coefficients of the M × M (M = 2d + 1) matrix equation are set (Eq. ( 1)). As done in the homogeneous and non-dissipative case, the discrete nodal pressure is expressed on the basis of

n p r = 2 2d + η + ζ 2d m=1 n I m r + η n I 2d+1 r , (17a)
n p r = 2d m=1 n T m r n I m r + n T 2d+1 r n I 2d+1 r = 2d+1 m=1 n T m r n I m r . (17b)
For this demonstration, the vertical velocity gradient and the atmospheric absorption loss are considered as time-independent. Thus, the derivation of the equivalent wave equation starts by writing Eq. (17a) at the time step n + 1:

2d + η + ζ 2 n+1 p r = 2d m=1 n+1 I m r + η × n+1 I 2d+1 r . ( 18
)
Applying twice the Eqs. ( 4) and ( 6) leads to the pressure scheme:

2d + η + ζ 2 n+1 p r + 2d + η -ζ 2 n-1 p r = d m=1 n P (j1+δm1,...,j d +δ md ) + n P (j1-δm1,...,j d -δ md ) + η n P r . (19)
Finally, using Taylor expansions (as before for the homogeneous case), the terms multiplied by ζ complete the scheme presented in section 3. They approximate the first-order derivative of the pressure and the equivalent wave equation solved by the TLM model in a dissipative inhomogeneous medium is then:

1 c 2 TLM ∂ 2 p ∂t 2 -∇ 2 p + ζ (j,l) 2∆ 2 2∆t ∂p ∂t = O ∆t 4 ∆ 2 + O ∆t 3 ∆ 2 + O ∆ 2 , ( 20
)
with c TLM (x r ) = 2 η r + 2d ∆ ∆t . (21)
The meteorological effects are therefore modelled by defining an effective velocity c eff given by [5]:

160 c eff = γRT + w.u, (22)
with γ the heat capacity ratio, R the gas constant, T the temperature in Kelvin, w the horizontal wind vector and u the unit vector of the direction of sound propagation. Since the condition c TLM = c eff must be fulfilled to model the acoustic wave equation, the parameter η is expressed as:

η r = 2∆ 2 ∆t 2 √ γRT + w.u 2 -2d. ( 23
)

Stability analysis

The introduction of a plane wave solution of the form

P = P 0 exp i d l=1 k l x l -ωt
in Eq. ( 20) allows to express the norm of the wave vector k = [k 1 , . . . , k d ] T as:

k 2 = ω 2 c 2 TLM 1 + i 2ζ (η + 2d)ω∆t ,
which can be reformulated using a first order Taylor expansion:

k ≈ k 0TLM + iα TLM + O(ζ 2 ), (24)
with k 0TLM = ω c TLM the norm of the wave vector in the non-dissipative case and

α TLM = ζ 2(η + 2d)∆ .
Out of these results, it can be observed that the solution of Eq. ( 20) is an attenuated wave because its amplitude decrease with the distance. However, the dependence of the absorption coefficient α TLM on the parameter η shows that there is a coupling between the heterogeneity model and the absorption. This effect is critical for the validity of the model because the atmospheric absorption effect modeled by the TLM method will be affected by the wind conditions, which should not be the case.

Furthermore, as mentioned in [START_REF] Tsuchiya | Numerical simulation of sound wave propagation with sound absorption in time domain, 13th[END_REF], the theoretical absorption coefficient α is frequency-dependent whereas α TLM is not. Hence, this modeling technique of the atmospheric absorption is not suitable for broadband sound prediction in outdoor environments. However, post-filtering the numerical results in order to artificially apply the atmospheric absorption is still possible to compensate this missing feature in the TLM model. If only the heterogeneity of the sound velocity is considered (i.e ζ = 0), the insertion of a discretized plane wave in Eq. ( 19) leads to the dispersion equation:

ω = 1 ∆t arccos      2 d m=1 cos(k m ∆ ) + η r η r + 2d      . ( 25
)
According to Eq. ( 25) the condition for the scheme to be stable is that the argument of the arccosinus ∈ [-1, 1], which is met only if η ≥ 0 because the solution η ≤ -2d is unacceptable. The presence of the η r coefficient in Eq. ( 25) also shows that the scheme is dispersive. In order to characterize this dispersion, it is necessary to estimate the range of values taken by η r .

For a given simulation with valid results until a maximal frequency f max , the spatial discretization should be taken as: ∆ =

λ min N = c min f max N
, with c min the minimal sound velocity in the propagation medium. Then, in order to satisfy the stability condition:

∆ √ dc eff ≥ ∆t, and ∆t = ∆ √ dc max , (26)
taking c max as the maximal sound velocity in the propagation medium to respect the stability condition in the worst case. As in Ref. [START_REF] Guillaume | Application of the transmission line matrix method for outdoor sound propagation modelling -Part 1: Model presentation and evaluation[END_REF], linear vertical gradients of temperature (∂T /∂z= ± 0, 35 C o .m -1 ) and wind speed (∂w/∂z= ± 0, 2 s -1 ) are assumed. Here the altitude z=0 is just an artifact to set T (z=0) = 20 C o and w(z= 0) = 0 m.s -1 , there is no ground effect considered in the stability analysis. More realistic profiles, such as linear-logarithmic (called generally "lin-log") vertical velocity profiles, could be used but, given the number of points considered, they would not make a significant difference in the results [5]. Then, the values of η are calculated for downward (i.e downwind with a positive temperature vertical gradient) and upward refraction conditions (i.e upwind with a negative temperature vertical gradient). These examples are the worst case scenarios since they present the widest ranges of values for c eff . These situations give that η ∈ [0; 4.82 × 10 -1 ] for downward atmospherical conditions and η ∈ [0; 5.22 × 10 -1 ] for upward atmospherical conditions.

Figs. 8 and9 show that the dispersion error is maximal when η is maximal (i.e c TLM = min (c eff (z)) , z ∈ [z min , z max ]). The main difference with the homogeneous case is that, for θ = π 4 + β π 2 (β ∈ Z), the scheme presents a dispersion error (≈ 0, 4% for N = 5). However, theses plots also highlight that this formulation of the TLM model for inhomogeneous media leads to a maximal dispersion error at the bottom part of the grid in the downward case and at the upper part of the grid in the upward case. For long range sound propagation simulation, it is important to keep in mind these limitations of the TLM model. Table 1 displays the minimum and maximum relative dispersion errors for the downward case, which correspond to the errors in the diagonal and the axial direction of the mesh, respectively. The values show that the local errors are acceptable for N=5 and negligible for N=10 and they are close to those obtained in the upward case. However, these are local errors around a node and it is necessary to study their impact at a larger scale in order to oberve the propagation of the dispersion error. Moreover, it is known that dispersion affects the group speed more than the phase speed [START_REF] Trefethen | Group Velocity in Finite Difference Schemes[END_REF], so the following section aims to analyze its effects on this quantity through a numerical long-range propagation experiment.

Numerical experiment

To analyze and quantify the effects of the dispersion highlighted theoretically in sections 3 and 4, a numerical experiment is carried out. Indeed, determining the limits of the TLM predictions at high frequencies according to the distance from a sound source is critical when considering long range sound propagation. The evaluation of these limits would enable the characterization of a critical distance (source-receiver) for a given number of points per wavelength at the maximal frequency. It will set the range at which the dispersion error on the results is too high to consider them relevant.

Simulation setup

As highlighted through Figs. 4, 8 and 9, a 2D-cartesian mesh implies a higher dispersion error on the axial direction than along the diagonal direction. To analyze this effect in terms of sound pressure level predictions, the setup illustrated in Fig. 10 is simulated. It is composed of a unique sound source surrounded by two arrays of eleven receivers located along the horizontal and diagonal (θ = π 4 ) directions of the mesh. Long-distance propagation becomes challenging when the frequency reaches high values. The specifications of a configuration used in the following results is given by Table 2 to get an overview of the computing ressources needed. To process the approximately 2700 million points to be solved and support the associated memory load, the TLM model was developed in OpenCL programming language to allow the parallelization of calculations on graphics processing units (GPUs), the parameterization of the simulations being managed on the basis of Python scripts [START_REF] Guillaume | Optimized transmission line matrix model implementation for graphics processing units computing in built-up environment[END_REF]. In the case of the example, the implementation avoids the storage of a 10 Gb matrix and allows to perform the simulation in around 3 hours using three Nvidia GeForce cards GTX 1080 Ti. ∆ Ny × Nz are the number of integration points in the respective axes. The memory usage is given for a float32 numpy array. tcomp is the computation time needed for this configuration.

[m] ∆t[s] ∆ mic [m] N y × N z [-] Memory [Gb] t comp [s] r max /λ min [-] 1.43e-
Note that a comparison of the computational costs of the TLM model with those of more common methods is not straightforward. Indeed, even though the TLM is equivalent to a 2-point FDTD scheme in space and time, its implementation is not based on an explicit scheme of the form P n+1 = f (P n ). The decomposition into 2 × d pulses at each point allows the calculations to be divided into parallelizable computational subdomains as described in [START_REF] Guillaume | Optimized transmission line matrix model implementation for graphics processing units computing in built-up environment[END_REF]. This method significantly reduces the computation time but multiplies by 2 × d the memory usage.

Sound source

The omnidirectional point source used in this experiment is a time-dependent Gaussian pulse, namely:

S(r s , n) = exp -π 2 (f c n∆t -1) 2 , ( 27
)
with a parameter f c = f max /2 which guarantees that the frequency support of the source reaches the maximal frequency f max .

Analytical solution

The simulation results are compared with the analytical solution in free-field for a monopole source [START_REF] Duan | The sound power output of a monopole source in a cylindrical pipe containing area discontinuities[END_REF], defined as a spherical source with a Gaussian mass flow [START_REF] Bruneau | Fundamentals of Acoustics[END_REF]: ∆ ĥ(r, ω) + k 2 w ĥ(r, ω) = -δ(r -r s ),

q(r, t) = 4π 2 v(r = r 0 + ) exp(-π 2 (f c t -1) 2 )δ(r -r 0 )H(t) = Q 0 (t)δ(r -r 0 ), (28)
with ĥ(r, ω) = P (r, ω)

-iωρ 0 Q0 (w) and k w = ± w c .
The known solutions of this equation in free-field are the following Green functions: in 3D.

Ĝd (r, ω) =      1 4i H (1) 0 (-k w r) in 2D, H (
The acoustic pressure field can finally be expressed in the frequency and time domains respectively by: P (r, ω) = -iωρ 0 Q0 (w) Ĝd (r, ω)

and p(r, t) =T F -1 -iωρ 0 Q0 (w) Ĝd (r, ω) .

This analytical solution is then considered as a reference to estimate the validity of the model and the errors among the numerically predicted pressure fields.

Results

This section aims at analyzing the results of the numerical experiment with a broadband sound source (f max = 2400 Hz) comparing the case N = 10 and N = 5 points per wavelength. In Fig. 11, the envelope pressure of the signal at the receivers is displayed according to the reduced time (t n -r c 0 ) which allow an observation of the flattening of the signal along the propagation distance. The effect of the dispersion along the horizontal array of microphones is directly noticeable: the results displayed on Fig. 11 show an apparent decrease of the group speed while the propagation distance increases. For a better visualisation of this delay, the vertical lines on Fig. 11 indicate the moment when 95 % of the signal energy has reached the microphone. Only the N = 10 case is shown below because a coarser grid brings longer delays but the overall behavior of the pressure envelopes is equivalent. This phenomenon is confirmed by the Figs. 12(a) and 12(c) which shows a maximal relative error of 0.7% and 2.9% compared to the analytical solution on the horizontal array of microphones. However, even if the group speed is impacted, the Figs. 12(b) and 12(d) shows that the calculated equivalent sound pressure levels (Leqs) are close to the analytical solution which confirms that the TLM is a dispersive numerical method but not a dissipative one. The simulated Leqs are therefore not affected in free-field because they are calculated with windowed integrals, which integrate the whole signal energy.

To better understand the frequency behavior of the dispersion, a one-third octave band decomposition has been performed on the microphones signals and the relative group speed error has been calculated for each band. The results of this post-processing are displayed on Figs. 13 and 14 and they imply that the high frequency contents are more delayed compared to the analytical solution and to the results from the diagonal array. The difference between the N = 5 and the N = 10 case is significant and shows that the impact of the number of points per wavelength on the error is prevailing on the propagation of the error along the source-receiver distance.

These plots highlight that the inherent dispersion of the TLM model is impacting free-field long distance simulation results. Even if the impact is below 2% compared to the analytical solution at 450 m and N = 10 for 2400 Hz, further investigation is needed to characterize the effect of dispersion when a ground reflection is simulated. Indeed, if the high frequency components of the signal are delayed, interfence patterns may occur between the direct field pulse and the reflected signal and consequently impact the energy of the signal received at the microphones. In this case, the TLM dispersion could impact the simulated sound pressure levels which is critical in the field of sound propagation predictions in the outdoor environment.

Conclusion

The usage of the TLM model as a solver for long-range outdoor sound propagation was discussed through both theoretical and numerical aspects. The general formulation of the method was renewed and a proper way of linking it to the wave equation was proposed. Thanks to these developments, the inherent 320 anisotropic dispersion of the model has been highlighted in two different cases and the non-validity of the usual atmospheric absorption has been pointed out. Finally, a numerical experiment has been set to characterize the effects of the dispersion on long-range free-field outdoor propagation results. From the results of this experiment and the theoretical study, the conclusion that the TLM results are valid in this configuration can be drawn provided that the number of point per wavelength is chosen wisely. Indeed, if N is too high, it will lead to significant computational costs and if N is too low, the dispersion error will reduce the frequency range of validity of the results compared to the maximal simulated frequency. However, the observed effects on the group speed of the test signals in the main mesh directions must be taken into account when using the TLM. In the case of predictions for sound propagation in outdoor environments, equivalent sound pressure levels are the primary quantity of interest for the evaluation of sound disturbances. For now, the importance of TLM dispersion in presence of boundary conditions such as natural ground is unknown. Thus, future work needs to be done regarding TLM for modeling realistic long-range sound fields.

Figure 1 :

 1 Figure 1: Phenomenological representation of an incident pulse reaching a node: (a) as a medium discontinuity, (b) on a 2D grid. The convention adopted for the lines numerotation is from left to right and from bottom to top.

  80

Figure 4 :

 4 Figure 4: Relative errors on the phase speed [%] related to the model anisotropy depending on the number N: (a) N=2, (b) N=5 and (c) N=10.

Figure 5 :

 5 Figure 5: Node transmission lines 2D environment example for an inhomogeneous and dissipative medium.

Figure 6 :

 6 Figure 6: Electrical analogy scheme along regular transmission lines in the case of an inhomogeneous dissipative medium.

Figure 7 :

 7 Figure 7: Electrical analogy scheme along the transmission line 2d + 1 in the case of an inhomogeneous dissipative medium.

  . The parameter ζ is calculated from the atmospheric absorption coefficient α [m -1 ]. It must be pointed out that, as the α coefficient, the ζ parameter is frequency-

Figure 8 :

 8 Figure 8: Modeling of downward refraction conditions for N = 5 : (a) evolution of η according to the height z; (b) relative errors on the phase speed [%] according to z and θ.

Figure 9 :

 9 Figure 9: Modeling of upward refraction conditions for N = 5 : (a) evolution of η according to the height z; (b) relative errors on the phase speed [%] according to z and θ.

Figure 10 :

 10 Figure 10: Dispersion analysis setup: the source in the center and microphones are positioned along the horizontal and diagonal (θ = π 4 ) directions.

  H(t) being the Heaviside function and δ(r -r 0 ) the spatial Dirac density such that Ω δ(r -r 0 )dΩ = 1. Applying the Fourier transform to the sound wave equation leads to the Helmoltz equation:

  1) 0 being the Hankel function of the first kind; exp(-ikwr) 4πr

Figure 11 :

 11 Figure 11: Normalized pressure envelopes at the microphones as dependent upon reduced time, for N=10 grid points per wavelength. Vertical lines indicates the arrival of 95 % of the signal energy.

Figure 12 :

 12 Figure 12: Comparison between the signals received at the microphones and the analytical solution depending on the propagation distance (a) & (c) magnitude of the errors on the group velocities for N = 10 and N = 5 respectively; (b) & (d) absolute errors on Leqs for N = 10 and N = 5 respectively.

Figure 13 :

 13 Figure 13: group speed error plots depending on the one-third octave band decomposed frequencies and the distance to the source for N = 10 for the horizontal array.

Figure 14 :

 14 Figure 14: group speed error plots depending on the one-third octave band decomposed frequencies and the distance to the source for N = 5 for the horizontal array.

Table 1 :

 1 Minimum and maximum values of the relative error on phase speed for various (η, N) configurations, downward case.

	η	N	2	5	10
		min [%]	max [%]	min [%] max [%]	min [%] max [%]
	4.89e-01	3.85e+00 3.03e+01 3.85e-01 3.87e+00 9.12e-02 9.26e-01
	3.16e-01	2.52e+00 2.86e+01 2.49e-01 3.60e+00 5.89e-02 8.61e-01
	1.53e-01	1.24e+00 2.70e+01 1.21e-01 3.35e+00 2.85e-02 8.01e-01
	0.00e+00	1.56e-04	2.55e+01 1.74e-05 3.12e+00 4.18e-06 7.44e-01

Table 2 :

 2 Configuration example for a setup with fmax = 2400 Hz, N = 10, rmax = 450 m.

	02 2.95e-05	45	56476 × 48746	10.22	11143	3176

Acknowledgments

The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data.


Sophie Mockly

Élisabeth Houbron

Hervé Seitz
email: herve.seitz@igh.cnrs.fr

A rationalized definition of general tumor suppressor microRNAs excludes miR-34a

Keywords: microRNA, miR-34, proliferation, tumor suppressor

While several microRNAs (miRNAs) have been proposed to act as tumor suppressors, a consensual definition of tumor suppressing miRNAs is still missing. Similarly to coding genes, we propose that tumor suppressor miRNAs must show evidence of genetic or epigenetic inactivation in cancers, and exhibit an anti-tumorigenic (e.g., anti-proliferative) activity under endogenous expression levels. Here we observe that this definition excludes the most extensively studied tumor suppressor candidate miRNA, miR-34a. In analyzable cancer types, miR-34a does not appear to be down-regulated in primary tumors relatively to normal adjacent tissues. Deletion of miR-34a is occasionally found in human cancers, but it does not seem to be driven by an anti-tumorigenic activity of the miRNA, since it is not observed upon smaller, miR-34a-specific alterations. Its anti-proliferative action was observed upon large, supra-physiological transfection of synthetic miR-34a in cultured cells, and our data indicates that endogenous miR-34a levels do not have such an effect. Our results therefore argue against a general tumor suppressive function for miR-34a, providing an explanation to the lack of efficiency of synthetic miR-34a administration against solid tumors.

Introduction

Tumor suppressors are genes whose activity antagonizes tumorigenesis. Consequently, they are frequently silenced, either by germline-inherited or somatic mutation, or otherwise inactivated, in cancers [START_REF] Green | Recessive mechanisms of malignancy[END_REF]. Mechanistically, tumor suppressors typically mediate cellular environment-induced inhibition of cell proliferation, therefore exhibiting anti-proliferative activity under their natural expression levels: a gene displaying cytotoxic or cytostatic activity only when inappropriately overexpressed is therefore excluded from that definition [START_REF] Weinberg | Tumor suppressor genes[END_REF]. miRNAs are small regulatory RNAs, guiding their effector proteins to specific target RNAs, which are repressed by various mechanisms (target RNA degradation and translational inhibition) [START_REF] Iwakawa | The functions of MicroRNAs: mRNA decay and translational repression[END_REF]. Targets are recognized by sequence complementarity, with most targets bearing a perfect match to the miRNA "seed" (nt 2-7) [START_REF] Bartel | MicroRNAs: target recognition and regulatory functions[END_REF]. Such a short binding motif makes miRNA/target binding poorly specific, and more than 60% of human genes are predicted to be targeted by at least one miRNA [START_REF] Friedman | Most mammalian mRNAs are conserved targets of microRNAs[END_REF].

Because such gene regulators can act in signal transduction cascades, they may participate in tumorsuppressive pathways. A consensual definition for "tumor suppressor miRNAs" is still lacking, with some tentative definitions being based on miRNA down-regulation in cancer cells [START_REF] Zhang | microRNAs as oncogenes and tumor suppressors[END_REF], on the targets' annotation [START_REF] Wong | DNA methylation of tumor suppressor miRNA genes: a lesson from the miR-34 family[END_REF], or both [START_REF] Adams | The tumor-suppressive and potential therapeutic functions of miR-34a in epithelial carcinomas[END_REF]. We rather propose to follow the initial definition of tumor suppressors [START_REF] Weinberg | Tumor suppressor genes[END_REF], considering that there is no reason to particularize miRNAs among other types of tumor suppressors.

We thus advocate for the following definition of tumor suppressor miRNAs: (i) there is evidence for their frequent inactivation in cancer (either by genetic or epigenetic alteration; potentially only in specific cancer types); and (ii) they inhibit tumorigenesis (e.g., by repressing cell proliferation) under their endogenous expression level, rather than upon unrealistic overexpression.

We applied this definition to interrogate the status of the most highly-studied tumor suppressor candidate miRNA, miR-34a. It is a member of the miR-34 family, comprising six members in human and in mouse: miR-34a, miR-34b, miR-34c, miR-449a, miR-449b and miR-449c (Supplementary Figure S1). The three miR-34a/b/c subfamily members are transcriptionally controlled by the p53 tumor suppressor, which suggested that these miRNAs could participate in the tumor suppressive activity of the p53 network [START_REF] He | A microRNA component of the p53 tumour suppressor network[END_REF][START_REF] Bommer | p53-mediated activation of miRNA34 candidate tumor-suppressor genes[END_REF][START_REF] Chang | Transactivation of miR-34a by p53 broadly influences gene expression and promotes apoptosis[END_REF][START_REF] Corney | MicroRNA-34b and MicroRNA-34c are targets of p53 and cooperate in control of cell proliferation and adhesion-independent growth[END_REF][START_REF] Tarasov | Differential regulation of microRNAs by p53 revealed by massively parallel sequencing: miR-34a is a p53 target that induces apoptosis and G1-arrest[END_REF][START_REF] Tazawa | Tumor-suppressive miR-34a induces senescence-like growth arrest through modulation of the E2F pathway in human colon cancer cells[END_REF][START_REF] Raver-Shapira | Transcriptional activation of miR-34a contributes to p53-mediated apoptosis[END_REF]. Indeed, the miR-34a member is down-regulated or lost in various cancer models (tumor samples or transformed cell lines) relatively to normal samples [START_REF] He | A microRNA component of the p53 tumour suppressor network[END_REF][START_REF] Chang | Transactivation of miR-34a by p53 broadly influences gene expression and promotes apoptosis[END_REF][START_REF] Welch | MicroRNA-34a functions as a potential tumor suppressor by inducing apoptosis in neuroblastoma cells[END_REF][START_REF] Bommer | p53-mediated activation of miRNA34 candidate tumor-suppressor genes[END_REF][START_REF] Tazawa | Tumor-suppressive miR-34a induces senescence-like growth arrest through modulation of the E2F pathway in human colon cancer cells[END_REF][START_REF] Lodygin | Inactivation of miR-34a by aberrant CpG methylation in multiple types of cancer[END_REF][START_REF] Gallardo | miR-34a as a prognostic marker of relapse in surgically resected non-small-cell lung cancer[END_REF][START_REF] Wiggins | Development of a lung cancer therapeutic based on the tumor suppressor microRNA-34[END_REF][START_REF] Corney | Frequent downregulation of miR-34 family in human ovarian cancers[END_REF]. This observation suggested that the inactivation of miR-34a is involved in tumorigenesis, and that other family members (miR-34b and c, miR-449a, b and c) could not compensate for this loss. miR-34a is therefore widely perceived as a general tumor suppressor, whose inactivation is involved in a variety of cancer types [START_REF] Slack | The Role of Non-coding RNAs in Oncology[END_REF]. Yet miR-34a -/-, miR-34b -/-, miR-34c -/-triple knock-out mice do not exhibit obvious defects in p53-dependent proliferation control or in tumor suppression [START_REF] Concepcion | Intact p53-dependent responses in miR-34-deficient mice[END_REF]. And, while pre-clinical studies in mice gave encouraging results (reviewed in [START_REF] Bader | miR-34 -a microRNA replacement therapy is headed to the clinic[END_REF][START_REF] Hong | Phase 1 study of MRX34, a liposomal miR-34a mimic, in patients with advanced solid tumours[END_REF]), administration of a synthetic miR-34a to human patients with solid tumors failed to repress tumor growth reproducibly [START_REF] Beg | Phase I study of MRX34, a liposomal miR-34a mimic, administered twice weekly in patients with advanced solid tumors[END_REF]. An alternate administration regimen (allowing increased drug exposure) did not clearly improve clinical outcomes, while triggering poorly-understood, severe adverse effects [START_REF] Hong | Phase 1 study of MRX34, a liposomal miR-34a mimic, in patients with advanced solid tumours[END_REF].

Materials and Methods

Analysis of miR-34a expression and integrity in human cancers

miRNA expression data was downloaded from the GDC portal on April 29, 2021. Cancer types where at least 10 cases were available (with Small RNA-Seq data from normal solid tissue and primary tumor for each case) were selected, and depth-normalized read counts were compared between normal tissue and tumor for each case. The heatmap shown on Figure 1A shows the median log-ratio between tumor and normal tissue, with non-significant changes (calculated with the Wilcoxon test, FDR-adjusted for multiple hypothesis testing) being colored in white. miRNA gene ploidy data was downloaded from the GDC portal on March 4, 2021. Erroneous miRNA gene coordinates were corrected using information from miRBase. For the heatmap shown on Figure 1B, the percentage of cases with miRNA gene loss (either homo-or heterozygous) was evaluated for each miRNA, selecting cancer types where ploidy was determined in at least 100 cases. miRNA sequence variation data was downloaded from the GDC portal on February 24, 2021. SNP location was intersected with miRNA hairpin and mature miRNA coordinates from miRBase (as well as with miRNA seed coordinates, defined as nt 2-7 of the mature miRNA). For the heatmaps shown on Supplementary Figure S2, the percentage of cases with sequence variations in miRNA genes (hairpin, mature or seed sequences) is displayed, selecting cancer types with at least 100 analyzed cases.

For each of these heatmaps, miRNAs and cancer types were clustered with the heatmap.2 command with the R software.

CRISPR/Cas9-mediated mutagenesis

Four sgRNAs were designed using CRISPOR (http://crispor.tefor.net/ [START_REF] Concordet | CRISPOR: intuitive guide selection for CRISPR/Cas9 genome editing experiments and screens[END_REF]) to target each side of the human pre-mir-34a sequence, and cloned into an expression plasmid for S. pyogenes Cas9 (pSpCas9(BB)-2A-GFP plasmid (PX458), a gift from Feng Zhang [START_REF] Ran | Genome engineering using the CRISPR-Cas9 system[END_REF]; Addgene plasmid #48138; http://n2t.net/addgene:48138; RRID:Addgene_48138). Targeting efficiency of each plasmid was estimated by Sanger sequencing of the targeted locus in transfected HCT-116 cells, and analyzed with the Synthego ICE Analysis online tool (https://ice.synthego.com/#/). Mutagenesis was performed using the most efficient sgRNA sequence on each side of the targeted locus (AAGCTCTTCT-GCGCCACGGTGGG and GCCGGTCCACGGCATCCGGAGGG; PAM sequences in bold; also see Supplementary Figure S5). HCT-116 (ATCC® cat. #CCL247) and HAP1 (Horizon Discovery cat. #C631) cells were grown till 80% confluency and transfected with the two plasmids (15 µg each) following the protocol for Lipofectamine 2000 Transfection Reagent (Thermo Fisher Scientific). After 24 hours, Cas9-GFPexpressing single cells were isolated in 96-well plates by flow cytometry on a BD FACSMelody (Becton Dickinson), then grown for 10 days. Homozygous wild-type and mutant clones were first tested by discriminative PCRs (with primer pairs ACTTCTAGGGCAGTATACTTGCT and GCT-GTGAGTGTTTCTTTGGC; and TCCTCCCCACATTTCCTTCT and GCAAACTTCTCCCAGC-CAAA), and eventually validated by Sanger sequencing of their miR-34a locus. For the HAP1 cell line, mutagenesis efficiency was so high that we were unable to isolate wild-type clones after cotransfection of sgRNA-carrying PX458 plasmids. Wild-type clones were therefore generated by transfection of HAP1 cells with a plasmid expressing SpCas9-HF1 variant but no sgRNA (the VP12 plasmid, a gift from Keith Joung [START_REF] Kleinstiver | High-fidelity CRISPR-Cas9 nucleases with no detectable genome-wide offtarget effects[END_REF]; Addgene plasmid #72247 ; http://n2t.net/addgene:72247; RRID:Addgene_72247), and went through the same isolation and selection process as mutant clones.

RNA extraction

Cells plated in 10 cm Petri dishes were lysed and scrapped in 6 mL ice-cold TRIzol™ Reagent (Invitrogen) added directly to the culture dish after removal of the growth medium, and mixed with 1.2 mL of water-saturated chloroform. Samples were homogenized by vigorous shaking for 1 min and centrifuged for 5 min at 12,000 g and 4°C to allow phase separation. The aqueous phase was transfered in a new tube and mixed with 3 mL isopropanol for precipitation. After a 10 min incubation at room temperature, samples were centrifuged for 10 min at 12,000 g and 4°C and the supernatant was removed. The RNA pellet was washed with 6 mL of 70% ethanol and samples were centrifuged for 5 min at 12,000 g and 4°C. After complete removal of ethanol, the RNA pellet was resuspended in 20 µL RNase-free water and the quantity of total RNA was determined by spectrophotometry on a NanoDrop ND-1000.

Small RNA-Seq

Total RNA of each cell line was extracted 48 h after seeding and quality was assessed on electrophoretic spectra from a Fragment Analyzer (Agilent), analyzed with the PROSize software (v. 3.0.1.6).

Libraries were prepared using NEXTflex™ Small RNA-Seq Kit v3 (Bioo Scientific) following the manufacturer's instructions. Libraries were verified by DNA quantification using Fragment Analyzer (kit High Sensitivity NGS), and by qPCR (ROCHE Light Cycler 480). Libraries were sequenced on Illumina NovaSeq 6000 using NovaSeq Reagent Kit (100 cycles). RNA quality assessment, library preparation, validation and sequencing were performed by the MGX sequencing facility.

Adapters ended with 4 randomized nucleotides in order to reduce ligation biases. Because of the sequencing design, the adapter sequence (5 ́GTTCAGAGTTCTACAGTCCGACGATCNNNN 3 ) appears at the beginning of the read sequence, and the final 4 nucleotides of the read are the initial randomized nucleotides of the other adapter, whose other nucleotides are not read. Hence small RNA reads can be extracted from the fastq files with the following command:

cutadapt -g GTTCAGAGTTCTACAGTCCGACGATCNNNN --discard-untrimmed -m 18 -M 30 \ $input_file.fastq | cutadapt -u -4 -

Cell transfection

Cells were transfected 24 hours after seeding either with a control duplex, siRNA against eGFP: 5 -GGCAAGCUGACCCUGAAGUdTdT-3 ́/ 5 -ACUUCAGGGUCAGCUUGCCdTdT-3 or with a hsa-miR-34a mimic duplex: 5 -P-UGGCAGUGUCUUAGCUGGUUGUU-3 ́/ 5 -P-CAAUCAGCAAGUAUACUGCCCUA-3 according to the protocol for Lipofectamine 2000 Transfection Reagent (Thermo Fisher Scientific).

Proliferation assays

Because the mere procedure of isolating and selecting mutated clones may artifactually select clones with exceptionally high proliferation rates, we applied the same isolation and selection procedure to wild-type clones, and we measured proliferation rates on several independent wild-type and mutant clones. Each cell line was seeded in 96-well plates (Figure 3C: in 4 replicates at 3 × 10 3 cells/well per time point; Figures 4A andB: in 6 replicates at 6×10 3 cells/well). From 24 hours after cell seeding or transfection, to 3 days later, the number of living cells was determined twice a day by CellTiter-Glo Luminescent Cell Viability Assay (Promega) according to the manufacturer's protocol and recorded with a TriStar LB 941 (Berthold Technologies). Linear regression of log-transformed cell counts relative to time and genotype (in R syntax: log-transformed cell counts ~time * genotype) or transfected duplex identity (log-transformed cell counts ~time * duplex identity) was used to measure doubling time and to estimate the significance of the effect of genotype or transfected duplex.

For Figure 3D andE, doxorubicin (Sigma-Aldrich) was diluted in molecular biology-grade water and 5-fluorouracil (5-FU) (Sigma-Aldrich) diluted in dimethyl sulfoxide (SigmaAldrich). In a preliminary experiment, half-maximal inhibitory concentration (IC50) was estimated after 72 h drug exposure: 7 × 10 -8 M and 8 × 10 -6 M for doxorubicin and 5-FU respectively. Cell lines were seeded in 3 replicates per drug concentration at 2.5 × 10 3 cells/well in 96-well plates. After 24 hours, culture medium was replaced with drug-containing medium (concentration range centered on the IC50 with 2.5× increments), or solvant-containing medium for untreated controls, and the number of living cells was determined 72 h later by CellTiterGlo Luminescent Cell Viability Assay (Promega). Cell counts were normalized to the mean cell number in untreated controls. Normalized cell number was fitted to an asymptotic model for each clone to assess the significance of the effect of genotype (using an analysis of variance to compare a model not informed by clone genotype, to a genotype-informed model).

Measurement of apoptosis

HCT-116 cells were seeded in 6-well plates in 3 replicates at 9 × 10 4 cells/well per condition. 72 hours after cell transfection, the number of apoptotic, dead and live cells was determined by FITC Annexin V/Dead Cell Apoptosis Kit with FITC annexin V and PI for flow cytometry (Invitrogen cat. #V13242), according to the manufacturer's protocol. Cells were analyzed by flow cytometry on a MACSQuant analyzer (Miltenyi) using the blue laser excitation (488 nm) with a 525/50BP filter. 10,000 singlet cells were measured per replicate, and apoptotic, dead and live populations were defined by FITC and PI thresholds pre-established with non-stained and mono-stained controls, and counted with the FlowJo Software (BD Biosciences).

miRNA quantification by RT-ddPCR

Reverse transcription of a specific miRNA in HCT-116 cells was performed on 10 ng total RNA using the TaqMan microRNA Reverse Transcription Kit (Thermo Fisher Scientific) in a total volume of 15 µL, according to the manufacturer's protocol, with miRNA-specific RT primers from the TaqMan MicroRNA Assay Kit (assay IDs for hsa-miR-34a-5p and miR-21b-5p are respectively 000426 and 000397). In order to ensure a precise ddPCR quantification, with similar numbers of positive and negative droplets in each sample, cDNA dilution factor was optimized for each experimental condition (cDNAs for miR-21 quantification were diluted 10×; cDNAs for miR-34 quantification were diluted:

100× for 1 nM-transfected samples, 1000× for 10 nM-transfected samples, undiluted for 0 nMtransfected samples and for samples shown in Figure 4D). ddPCR amplification of the cDNA was performed on 1.33 µL of each cDNA dilution combined with 1 µL of miRNA-specific 20X TaqMan MicroRNA Reagent containing probes and primers for amplification from the TaqMan MicroRNA Assay Kit (Thermo Fisher Scientific), 10 µL of 2X ddPCR Supermix for probes (no dUTP) (Bio-Rad), and 7.67 µL of molecular biology-grade water. Droplets were generated, thermal cycled and detected by the QX200 Droplet Digital PCR System (Bio-Rad) according to the ddPCR Supermix protocol and manufacturer's instructions. Data were extracted using QuantaSoft Pro Software (Bio-Rad).

Statistical analyses of RT-ddPCR data

miR-34a and miR-21 quantification was performed in 3 independent replicates, and cDNA counts were converted into numbers of miRNA molecules per ng RNA, considering dilution factors at each step in the RT-ddPCR process. Significance of the effect of transfected dose (for Figure 4C), doxorubicin treatment (for Figure 4D) and time (for both panels) was assessed by two-way ANOVA without interaction. Post-hoc pairwise t-tests were performed whenever the ANOVA found a significant effect for an explanatory variable.

Results

No evidence for miR-34a loss or inactivation in cancers

It is now possible to compare miRNA levels between tumors and normal adjacent tissues on a large collection of human cases which passed stringent, standardized quality controls [START_REF] Zhang | Uniform genomic data analysis in the NCI Genomic Data Commons[END_REF], allowing a rigorous assessment of miR-34a expression in tumorigenesis. Selecting every cancer type where miRNA expression is available for primary tumor and normal adjacent tissue, in at least 10 studied cases (n=20 cancer types), we did not find any cancer type where miR-34a was significantly down-regulated (Fig- ure 1A). Hence in this collection of cancer types, human primary tumors do not tend to under-express miR-34a, contradicting the notion that genetic or epigenetic silencing of miR-34a could participate in tumorigenesis.

Accordingly, genetic alterations affecting miR-34a are very rare in cancer: focusing on every cancer type for which gene-level copy number was measured in at least 100 cases (n=29 cancer types), we did not observe any tendency for the loss of miR-34a relatively to other miRNA genes (see Figure 1B). Similarly, we did not find any evidence for the selective mutation of the pre-miR-34a hairpin precursor sequence, mature miR-34a or the miR-34a seed in cancers (n=30 analyzed cancer types; Supplementary Figure S2). In contrast to miR-34a, 105 miRNA loci tend to be frequently lost in 19 cancer types (red area at the top left corner of the heatmap in Figure 1B; listed in Supplementary Table S1): these miRNAs are more convincing tumor suppressor candidates than miR-34a in this respect.

It could be argued that miR-34a inactivation by itself is insufficient to contribute to tumorigenesis, while it may play a role in a sensitized context, where additional, cooperative mutations may reveal the oncogenicity of miR-34a down-regulation. In that case, miR-34a inactivation could be enriched in just a subset of highly mutated cancers, and it would not be visible in the analyses shown in Figure 1 and Supplementary Figure S2. Yet, stratifying cases by cancer grade, we did not observe any tendency for the most aggressive tumors to inactivate miR-34a (Supplementary Figure S3), indicating that even the most sensitized tumors do not show any evidence of miR-34a inactivation.

Similarly, it is conceivable that miR-34a plays a tumor suppressive role only in the presence of functional p53, and the frequent mutation of p53 in the samples analyzed in Figure 1 may have obscured its behaviour in p53 +/+ tumors. But the selective analysis of cancer cases without any mutation in p53 gives a very similar result, without miR-34a being down-regulated in any analyzed cancer type (see Figure 2).

Hence the loss or mutation of miR-34a does not appear to be enriched in cancer. We note that miR-34a is located on cytogenetic band 1p36, which is often altered in a wide variety of cancers.

But our analyses suggest that the inactivation of miR-34a is not the actual driver for deletion selection -and because a convincing tumor suppressor is already known at 1p36 (the CHD5 gene [START_REF] Bagchi | CHD5 is a tumor suppressor at human 1p36[END_REF]), we propose that the occasional deletion of miR-34a in cancer is rather a consequence of its genomic proximity with such a real tumor suppressor. Accordingly, whenever a limited region of consistent deletion could be mapped in 1p36, that region excludes miR-34a (with the only exception of myelodysplastic syndromes, but with low experimental support): see Supplementary Figure S4.

The reported anti-proliferative action of miR-34a is artifactual

miR-34a has also been considered a tumor suppressor candidate on the basis of the apparent antiproliferative activity of miR-34 family miRNAs. Numerous studies in cultured cell lines indeed showed that miR-34 transfection inhibits cell proliferation, either by slowing down cell division or by increasing cell death [START_REF] Welch | MicroRNA-34a functions as a potential tumor suppressor by inducing apoptosis in neuroblastoma cells[END_REF][START_REF] He | A microRNA component of the p53 tumour suppressor network[END_REF][START_REF] Chang | Transactivation of miR-34a by p53 broadly influences gene expression and promotes apoptosis[END_REF][START_REF] Corney | MicroRNA-34b and MicroRNA-34c are targets of p53 and cooperate in control of cell proliferation and adhesion-independent growth[END_REF][START_REF] Tarasov | Differential regulation of microRNAs by p53 revealed by massively parallel sequencing: miR-34a is a p53 target that induces apoptosis and G1-arrest[END_REF][START_REF] Tazawa | Tumor-suppressive miR-34a induces senescence-like growth arrest through modulation of the E2F pathway in human colon cancer cells[END_REF][START_REF] Raver-Shapira | Transcriptional activation of miR-34a contributes to p53-mediated apoptosis[END_REF]. But miRNA over-expression generates false positives, raising the possibility that this reported anti-proliferative role is artifactual [START_REF] Mockly | Inconsistencies and Limitations of Current MicroRNA Target Identification Methods[END_REF]. We thus deleted the miR-34a gene in HCT-116 cells, where it has been proposed to be anti-proliferative by several independent studies [START_REF] He | A microRNA component of the p53 tumour suppressor network[END_REF][START_REF] Chang | Transactivation of miR-34a by p53 broadly influences gene expression and promotes apoptosis[END_REF][START_REF] Tazawa | Tumor-suppressive miR-34a induces senescence-like growth arrest through modulation of the E2F pathway in human colon cancer cells[END_REF] (mutagenesis strategy in Supplementary Figure S5). Deletion of the miR-34a locus eliminated 94% of the expression of the whole miR-34 family (Figure 3A andB). Our results do not show any significant difference in the growth rate of miR-34a -/-and wild-type clones (Figure 3C). We also prepared miR-34a -clones from the human haploid HAP1 cell line, where miR-34a is also not anti-proliferative (it is even slightly pro-proliferative; Supplementary Figure S6). It could be argued that miR-34a does not inhibit cell proliferation in unstressed conditions, while being anti-proliferative upon genotoxic stress. But we also failed to observe significant differences between wild-type and mutant clones under doxorubicin or 5-fluoro-uracil treatment (Figure 3D andE).

In agreement with published data, we did observe a strong reduction in cell proliferation when we transfected HCT-116 cells with large amounts (10 nM) synthetic miR-34a duplex (Figure 4A), but that effect was lost when transfecting 1 nM duplex (Figure 4B). Absolute miRNA quantification by RT-ddPCR shows that a 10 nM transfection over-expresses miR-34a by >8,000-fold in HCT-116 cells (and a 1 nM transfection over-expresses it by >490-fold), clearly demonstrating that such an experiment results in supra-physiological miRNA concentrations (Figure 4C). For comparison, we measured the increase in miR-34a expression in response to DNA damage: a 72 h treatment with doxorubicin at its IC50 concentration (7 × 10 -8 M in HCT-116 cells; Supplementary Figure S7) over-expresses miR-34a by only 4.7-fold (Figure 4D).

It could be argued that low doses of transfected miR-34a could induce apoptosis, and our measurements of cell viability may have missed early apoptotic cells, therefore under-estimating the cytotoxicity of low dose miR-34a. This interpretation is ruled out by the measurement of Annexin Vlabeled cells after transfection of various doses of miR-34a: physiological (picomolar range) doses do not appear to affect cell viability, and only the highest doses (>10 nM) lead to a visible decrease in cell viability (both through apoptosis and apoptosis-independent pathways): see Supplementary Figure S8.

Discussion

Our results show that the miR-34a gene is rarely inactivated in cancers, whether by deletion, mutation or any other kind of process affecting miRNA expression. Its occasional loss in some cancers is most likely due to its genomic proximity with an actual tumor suppressor, and miR-34a-specific mutations are not enriched in any cancer type with data available in the largest cancer genomic dataset available.

We also observed that the widely-assumed anti-proliferative role of miR-34a appears to be due to artifactual over-expression in cultured cells.

Of note, some authors have previously characterized the proliferative effect of miR-34 using genetic ablation rather than over-expression. In one study, mouse embryonic fibroblasts (MEFs) devoid of miR-34a/b/c appear to grow at the same rate than wild-type MEFs, except, transiently, for one early time-point [START_REF] Concepcion | Intact p53-dependent responses in miR-34-deficient mice[END_REF]. Mutation of miR-34a alone also appeared not to affect MEF proliferation [START_REF] Choi | miR-34 miRNAs provide a barrier for somatic cell reprogramming[END_REF].

In another study, genetic inactivation of the miR-34a gene in HCT-116 is reported to accelerate cell proliferation, in stark contrast with our own findings [START_REF] Navarro | miR-34 and p53: New Insights into a Complex Functional Relationship[END_REF]. Such discrepancy would deserve to be investigated, but unfortunately that published mutant cell line has been lost and it is no longer available from the authors (Dr. J. Lieberman, personal communication).

While the miR-34 family is believed to exert a tumor suppressive action in a diversity of cancers [START_REF] Slack | The Role of Non-coding RNAs in Oncology[END_REF], we observed that it is hardly expressed in cultured cell lines, primary tissues and body fluids (Supplementary Figures S9-S11). It could be argued that a low level of miR-34 expression is expected in normal tissues, where p53 is mostly inactive. But p53 is clearly not the only regulator for miR-34, and the expression of miR-34 does not mirror p53 activity [START_REF] Concepcion | Intact p53-dependent responses in miR-34-deficient mice[END_REF]. Current RNA detection technologies can be extremely sensitive, and they can detect miRNAs which are too poorly abundant to induce any clear change in target expression [START_REF] Mullokandov | High-throughput assessment of microRNA activity and function using microRNA sensor and decoy libraries[END_REF]. Hence we anticipate that in all the cell lines for which we analyzed miRNA abundance, and in most cells in the analyzed tissues, miR-34 family miRNAs are actually non-functional.

Yet we do not question the overall functionality of miR-34 miRNAs in vivo. Because that family is deeply conserved in evolution (shared between, e.g., vertebrates and insects), it certainly plays important biological functions, perhaps only in a small number of cells, or at very specific developmental stages, where its abundance would be high enough. In mouse, the miR-34 family is particularly expressed in lungs and testes [START_REF] Concepcion | Intact p53-dependent responses in miR-34-deficient mice[END_REF][START_REF] Song | miR-34/449 miRNAs are required for motile ciliogenesis by repressing cp110[END_REF]. Mutation of all 6 members of the miR-34 family causes severe ciliogenesis defects, leading to respiratory distress and impaired gametogenesis -translating into sterility and premature mortality [START_REF] Song | miR-34/449 miRNAs are required for motile ciliogenesis by repressing cp110[END_REF]. Unsurprisingly then, the most obvious biological functions of that miRNA family seem to take place in the tissues where miR-34 miRNAs are highly expressed, in contrast with the widely-accepted notion of their broad anti-tumorigenic activity.

While the original definition for tumor suppressors had been formulated with coding genes in mind, we consider that there is no objective reason for adopting a different definition for tumor suppressor miRNAs. In this view, the most heavily studied candidate tumor suppressor miRNA, miR-34a, does not appear to be a tumor suppressor. It remains formally possible that miR-34a inactivation is frequent in specific cancer types, distinct from those we could analyze in Figures 1 and2 and Supplementary Figure S2-3. In that case, miR-34a may be a tumor suppressor in these particular cancers, but rigorous investigation -while avoiding the pitfalls described above -of its impact on cell proliferation and tumorigenesis would be necessary to conclude so.

We confirmed that a large artificial over-expression (10 nM) of miR-34a indeed represses cell proliferation. It could be argued that this cytotoxic effect could provide the ground for an efficient anti-cancer treatment, no matter how un-natural it is. But the whole purpose of using natural tumor suppressors (e.g., miRNAs) is that they are expected to be well tolerated, because they already exist endogenously. Administering large amounts of cytotoxic agents to patients may indeed kill cancer cells -but it will also likely trigger unwanted adverse effects. In this view, synthetic miR-34a behaves similarly to existing anti-cancer drugs, which are based on exogenous molecules. It is therefore not surprising to observe a variety of adverse secondary effects when the MRX34 miR-34a mimic is administered to patients [START_REF] Beg | Phase I study of MRX34, a liposomal miR-34a mimic, administered twice weekly in patients with advanced solid tumors[END_REF][START_REF] Hong | Phase 1 study of MRX34, a liposomal miR-34a mimic, in patients with advanced solid tumours[END_REF]. More inocuous miRNA-based treatments may be possible, but they would have to rely on rigorously established tumor-suppressive activity of the endogenous miRNA.

Figure 1: mir-34a is not generally down-regulated or lost in cancers. A. miRNA abundance (normalized by the number of mapped miRNA reads) was compared between primary tumors and normal adjacent tissues. Only cancer types for which at least 10 cases were analyzed have been considered (n=20 cancer types; rows), and miRNAs with a null variance across cancer types were excluded (remaining: n=545 miRNAs; columns). For each miRNA/cancer type pair, the heatmap shows its median fold-change across all cases, with non-significant changes (FDR 0.05) being shown in white. log(fold-changes) larger than +8 or smaller than -8 were set to +8 or -8 respectively, for graphical clarity. B. Only cancer types for which at least 100 cases were analyzed have been considered (n=29 cancer types; rows), and miRNA genes whose ploidy could not be assessed were excluded (remaining: n=1,686 miRNAs; columns). For each miRNA/cancer type pair, the heatmap shows the percentage of cases with monoallelic or biallelic loss of the miRNA gene. Both panels: the column showing miR-34a data is magnified on the right margin (framed in black). "NOS": not otherwise specified.

Figure 2: No evidence for mir-34a inactivation in tumors with an intact p53 gene. Cancer samples analyzed in Figure 1A were stratified by the mutation status of the p53 gene, and only cases without any detected mutation in p53 were selected here (also selecting cancer types with at least 10 cases after this selection). Same conventions than in Figure 1A. miRNA abundance (normalized by the number of mapped miRNA reads) was compared between primary tumors and normal adjacent tissues. The column showing miR-34a data is magnified on the right margin (framed in black). log(fold-changes) larger than +5 or smaller than -5 were set to +5 or -5 respectively, for graphical clarity. "NOS": not otherwise specified. Linear modeling of log-transformed cell counts relative to time was used to measure doubling time (T d ), and to estimate the significance of the effect of genotype (p-value is given in the inset). Shaded areas represent the 95% confidence interval for theoretical future measurements. D, E. Cell number after 3 days of culture in presence of varying doses of (D) doxorubicin or (E) 5-fluoro-uracil (4 clones of each genotype were analyzed; 3 biological replicates for each drug concentration; mean +/-st. error is shown). Cell number was normalized to cell number count in untreated replicates. Normalized cell number was fitted to an asymptotic model for each clone (fitted models are represented by curves). In order to assess the significance of the effect of genotype, a naïve (non-informed by clone genotype) and a genotype-informed model were compared by an analysis of variance (p-value is indicated in the inset).

Figure 4: Supra-physiological transfection of miR-34a inhibits cell proliferation. Wild-type HCT-116 cells were transfected with 10 nM (panel A) or 1 nM (panel B) duplex (either a control siRNA duplex, or miR-34a/miR-34a* duplex) and grown in sub-confluent conditions. Means and standard errors of 6 biological replicates are represented by dots and error bars. Linear modeling of log-transformed cell counts relative to time was used to measure doubling time (T d ), and to estimate the significance of the effect of duplex identity (p-values are given in the inset; asterisks denote pvalue < 0.05, "n.s." indicates larger p-values). Shaded areas represent the 95% confidence interval for theoretical future measurements. C. Cellular abundance of miR-34a (red bars) or a control miRNA (miR-21; gray bars) 1 or 24 h after transfection of HCT-116 cells with 0, 1 or 10 nM miR-34a/miR-34a* duplex. D. HCT-116 cells were treated for 24 or 72 h with 7 × 10 -8 M doxorubicin, and their intracellular miR-34a and miR-21 were quantified by RT-ddPCR. Two-way ANOVA analysis shows that doxorubicin treatment has an effect on miR-34a levels (p=0.0013), and post-hoc pairwise t-tests find the effect significant only after 72 h exposure to the drug (p=0.0521 for 24 h exposure, p=0.00138 for 72 h exposure, indicated by "n.s." and "**" respectively). A similar two-way ANOVA analysis does not detect a significant effect of doxorubicin treatment on miR-21 levels (p=0.768). Panels C and D: Means and standard deviations of 3 biological replicates are represented by dots and error bars, respectively.

Figure 3 :

 3 Figure 3: miR-34 is not a general repressor of cell proliferation. A. miRNA quantification by Small RNA-Seq in a representative wild-type HCT-116 clone (x axis) and a representative miR-34a -/-clone (y axis). Right panel: magnification of the left panel. B. Cumulated abundance of miR-34 family members in the two clones. miRNAs are sorted vertically according to their abundance in the wild-type clone. C. Four wild-type and four miR-34a mutant clones were grown in sub-confluent conditions. Means and standard errors of 4 biological replicates are represented by dots and error bars.Linear modeling of log-transformed cell counts relative to time was used to measure doubling time (T d ), and to estimate the significance of the effect of genotype (p-value is given in the inset). Shaded areas represent the 95% confidence interval for theoretical future measurements. D, E. Cell number after 3 days of culture in presence of varying doses of (D) doxorubicin or (E) 5-fluoro-uracil (4 clones of each genotype were analyzed; 3 biological replicates for each drug concentration; mean +/-st. error is shown). Cell number was normalized to cell number count in untreated replicates. Normalized cell number was fitted to an asymptotic model for each clone (fitted models are represented by curves). In order to assess the significance of the effect of genotype, a naïve (non-informed by clone genotype) and a genotype-informed model were compared by an analysis of variance (p-value is indicated in the inset).

Acknowledgments

We thank G. Canal, C. Theillet, P.D. Zamore and members of the Seitz lab for critical reading of the manuscript, and A. Pélisson and K. Mochizuki for useful discussions. This research was supported by CNRS (salaries for É.H. and H.S., and running laboratory costs), by the university of Montpellier (PhD fellowship to S.M.), Cancéropôle GSO "Émergence" grant, and by Fondation ARC (grant "Projet" #PJA 20191209613 to the Seitz lab and 4 th year PhD fellowship #ARCDOC42019120001002 to S.M.). We thank the MGX facility (Biocampus Montpellier, CNRS, INSERM, Univ. Montpellier, Montpellier, France) for sequencing the Small RNA-Seq libraries.

Data and script availability

Deep-sequencing data has been deposited at SRA and linked to BioProject number PRJNA695193. Scripts, raw, intermediate and final data files are available at https://github.com/HKeyHKey/ Mockly_et_al_2022 and at https://www.igh.cnrs.fr/en/research/departments/genetics-development/ systemic-impact-of-small-regulatory-rnas#programmes-informatiques/. In particular, flow cytometry raw data has been deposited at https://github.com/HKeyHKey/Mockly_et_al_2022/ tree/main/Suppl_Figure_8/Flow_cytometry_raw_data.

Author contributions

S.M. and É.H. performed experiments; S.M. and H.S. performed computational analyses; S.M. and H.S. wrote the manuscript and prepared figures.

Conflict of interest

The authors do not declare any conflict of interest.