Arthur Marmin
email: arthur.marmin@univ-amu.fr

José Henrique De Morais Goulart

Cédric Févotte
email: cedric.fevotte@irit.fr

Joint Majorization-Minimization for Nonnegative Matrix Factorization with the β-divergence

Keywords: Nonnegative matrix multiplication (NMF), beta-divergence, joint optimization, majorization-minimization method (MM)

This article proposes new multiplicative updates for nonnegative matrix factorization (NMF) with the β-divergence objective function. Our new updates are derived from a joint majorization-minimization (MM) scheme, in which an auxiliary function (a tight upper bound of the objective function) is built for the two factors jointly and minimized at each iteration. This is in contrast with the classic approach in which a majorizer is derived for each factor separately.

Like that classic approach, our joint MM algorithm also results in multiplicative updates that are simple to implement. They however yield a significant drop of computation time (for equally good solutions), in particular for some β-divergences of important applicative interest, such as the quadratic loss and the Kullback-Leibler or Itakura-Saito divergences. We report experimental results using diverse datasets: face images, an audio spectrogram, hyperspectral data and song play counts. Depending on the value of β and on the dataset, our joint MM approach can yield CPU time reductions from about 13% to 86%

Introduction

Nonnegative matrix factorization (NMF) aims at factorizing a matrix with nonnegative entries into the product of two nonnegative matrices. It has found many applications in various domains which include feature extraction in image processing and text mining [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], audio source separation [START_REF] Smaragdis | Static and dynamic source separation using nonnegative factorizations: A unified view[END_REF], blind unmixing in hyperspectral imaging [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF][START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF], and user recommendation [START_REF] Hu | Collaborative filtering for implicit feed-back datasets[END_REF]. See [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation[END_REF][START_REF] Fu | Nonnegative matrix factorization for signal and data analytics: Identifiability, algorithms, and applications[END_REF][START_REF] Gillis | Nonnegative Matrix Factorization[END_REF] for overview papers and books about NMF.

Each application gives different interpretations to the factor matrices but the first factor is often considered as a dictionary of recurring patterns while the second one describes how the data samples are expanded onto the dictionary (activation matrix). The nonnegativity constraint only allows additive combination of the dictionary elements, yielding meaningful additive and sparse representations of the data.

Computing an NMF generally consists in minimizing a well-chosen objective function under nonnegativity constraints. A popular choice of objective function is the β-divergence, which is a continuous family of measures of fit parameterized by a single parameter β that encompasses the Kullback-Leibler (KL) or Itakura-Saito (IS) divergences as well as the common squared Euclidean distance [START_REF] Cichocki | Generalized Alpha-Beta divergences and their application to robust nonnegative matrix factorization[END_REF].

In the latter cases, the β-divergence is a log-likelihood in disguise for Poisson, multiplicative Gamma and additive Gaussian noise models, respectively.

The classic approach to NMF, and to NMF with the β-divergence in particular, consists in optimizing the two factors alternately, i.e., using two-block coordinate descent. Each of the two factors is then updated using Majorization-Minimization (MM), as described in [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF][START_REF] Kompass | A generalized divergence measure for nonnegative matrix factorization[END_REF][START_REF] Nakano | Convergence-guaranteed multiplicative algorithms for nonnegative matrix factorization with β-divergence[END_REF][START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF][START_REF] Yang | Unified development of multiplicative algorithms for linear and quadratic nonnegative matrix factorization[END_REF]: given one of the factors, a tight upper bound of the objective function is constructed and minimized with respect to (w.r.t) the other factor. This results in multiplicative updates that are simple to implement (with no hyperparameter to tune), have linear complexity per iteration, and that automatically preserve nonnegativity given positive initializations. By construction, MM ensures monotonicity of the objective function (non-increasing values), see [START_REF] Lange | MM Optimization Algorithms[END_REF][START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF] for tutorials about MM.

Thanks to its simplicity, the block-descent approach is dominant in matrix factorization and dictionary learning (using sometimes other block partitions such as columns or rows [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation[END_REF][START_REF] Gillis | Nonnegative Matrix Factorization[END_REF]) and very few works have addressed joint ("allat-once") optimization of the factors, though the latter approach may be more efficient. A notable exception in dictionary learning (real-valued factors, sparse activations, quadratic loss) is [START_REF] Rakotomamonjy | Direct optimization of the dictionary learning problem[END_REF]. In this work, the author employs an elegant non-convex proximal splitting strategy and shows that the joint approach is significantly faster than alternating methods without altering the quality of the obtained solution. In a similar spirit to [START_REF] Rakotomamonjy | Direct optimization of the dictionary learning problem[END_REF], [START_REF] Mukkamala | Beyond alternating updates for matrix factorization with inertial Bregman proximal gradient algorithms[END_REF] leverages the theoretical framework of [START_REF] Bolte | First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems[END_REF] to address matrix factorization (including NMF) with non-alternating updates. Their work relies on the generalization of Lipschitzcontinuity of the gradient (which does not hold jointly for both factors) to adaptive smoothness [START_REF] Bolte | First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems[END_REF]. Yet, their results only apply to the quadratic loss and Newton-like acceleration is crucial to obtain competitive results. Using tools from dynamical systems, the authors in [START_REF] Panageas | Convergence to second-order stationarity for non-negative matrix factorization: Provably and concurrently[END_REF] have derived a novel form of multiplicative updates which can run concurrently for each factor matrix at a given iteration. They have the additional benefit of ensuring the convergence to a local minimum instead of a mere critical point. However, their results are again limited to the quadratic loss. This also applies to [START_REF] Marumo | Majorization-minimizationbased Levenberg-Marquardt method for constrained nonlinear least squares[END_REF] where a Levenberg-Marquardt joint optimization method is described for NMF with the quadratic loss. In [START_REF] Vandecappelle | A second-order method for fitting the canonical polyadic decomposition with nonleast-squares cost[END_REF], the authors propose a joint second-order Newton-like algorithm for nonnegative canonical tensor decomposition with the β-divergence, which takes NMF as a special case. Their approach relies on approximations of the Hessian matrix, sometimes based on heuristics, and fails to provide an algorithm that universally works for every value of β.

Inspired by these works, we propose a joint MM approach to β-NMF and compare it to the classic block MM strategy. Our joint MM relies on a tight majorization of the objective function with respect to all its variables instead of using blocks of fixed variables. Iterative minimization of this joint upper bound results in new multiplicative updates that are simple but potentially more efficient variants of the classical multiplicative updates. This is particularly true when considering the quadratic loss and the KL or IS divergences for which further simplifications occur. We show in these cases that our update rules decrease the computation cost per iteration. It turns out that our joint upper bound coincides with the one derived in [START_REF] Takahashi | A unified global convergence analysis of multiplicative update rules for nonnegative matrix factorization[END_REF]. The latter is however employed for a different purpose, namely the convergence analysis of classic block MM for NMF, and not to design a new algorithm like we do (more details will be given in Section 3.3).

Our methodological results are supported by extensive simulations using datasets with different sizes arising from various applications in which NMF had a significant impact (face images, audio spectrogram, hyperspectral images, song play-counts). In most scenarii, the proposed joint MM approach leads to a reduction of the computing time ranging from 13% to 86% without incurring any loss in the precision of the approximation.

The article is organized as follows: Section 2 states the NMF optimization problem and summarizes the classic block MM approach. Section 3 first presents our proposed joint MM method for NMF and derives the new multiplicative updates. It then compares the joint and the block MM methods in term of computational complexity, and discusses the benefit of the joint approach. Comparative numerical simulations are presented in Section 4 and validate the efficiency of our approach. Finally, Section 5 draws conclusions.

Notation. R + is the set of nonnegative real numbers, and 1, N is the set of integers from 1 to N . Bold upper case letters denote matrices, bold lower case letters denote vectors, and lower case letters denote scalars. The notation [M] ij and m ij both stand for the element of M located at the i th row and the j th column. For a matrix M, the notation M ≥ 0 denotes entry-wise nonnegativity.

Preliminaries

Nonnegative matrix factorization

We aim at factorizing a F × N nonnegative matrix V into the product WH of two nonnegative factor matrices of sizes F × K and K × N , respectively. The rank value K is often chosen such that F K + KN ≪ F N , leading to a low-rank approximation of V. Given a desired value for the rank K, the factor matrices are obtained by solving the following optimization problem min

W,H≥0 D β (V | WH) , (1)
where D β is a separable objective function defined by

D β (V | WH) = F f =1 N n=1 d β (v f n | [WH] f n ) . (2)
Our measure of fit d β is the β-divergence [START_REF] Cichocki | Generalized Alpha-Beta divergences and their application to robust nonnegative matrix factorization[END_REF] given by

d β (x | y) =            x log x y -x + y if β = 1 x y -log x y -1 if β = 0 x β β(β-1) + y β β -xy β-1 β-1
otherwise.

(

) 3
The value of β is chosen according to the application context and the noise assumptions on V [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF]. The IS divergence, KL divergence and quadratic loss are obtained for β = 0, 1, 2, respectively.

Classic multiplicative updates

The classic method to solve Problem (1) consists in a two-block coordinate descent approach where each block is handled with MM. Namely, it alternately minimizes (W, H) → D β (V|WH) in W and in H using MM. The MM method is a two-step iterative optimization scheme [START_REF] Lange | MM Optimization Algorithms[END_REF][START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF]. At each iteration, the first step consists in building a local auxiliary function G that is minimized in the second step. The auxiliary function has to be a tight majorizer of the original objective function φ : E → R at the current iterate x, where E is the domain of G and φ. More precisely, it has to satisfy the following two properties:

(∀x ∈ E) G(x | x) ≥ φ(x) G(x | x) = φ(x) .
These properties ensure that any iterate x that decreases the value of G also decreases the value of φ. Indeed, for a given x, if we find x such that G(x|x) ≤ G(x|x), then the tight majorization properties induce the following descent lemma

φ(x) ≤ G(x | x) ≤ G(x | x) = φ(x) . (4)
Note that even if G is not minimized but only decreased in value, the descent property still holds.

The previous MM scheme can be applied alternately to the minimization of the two functions W → D β (V|WH) and H → D β (V|WH). These two functions are the sum of concave, convex, and constant terms. A convex auxiliary function can then be easily derived by using Jensen's inequality for the convex term and the tangent inequality for the concave term, see [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF] and the next section. This yields the following multiplicative updates

W ←-W.   (WH) .(β-2) .V H ⊤ (WH) .(β-1) H ⊤   .γ(β) H ←-H.   W ⊤ (WH) .(β-2) .V W ⊤ (WH) .(β-1)   .γ(β) , (5)
where . and / are the entry-wise multiplication and division, respectively, and γ(β) is a scalar defined as

γ(β) =            1 2-β if β ∈] -∞, 1[ 1 if β ∈ [1, 2] 1 β-1 if β ∈]2, +∞[ . (6)
Note that by construction, the matrices W and H resulting from the update [START_REF] Hu | Collaborative filtering for implicit feed-back datasets[END_REF] contain only positive coefficients if the input matrices V, W, and H are all positive. The nonnegativity of the iterates is thus preserved by the multiplicative structure of the updates given positive initializations.

Note that strict MM dictates that the individual updates of W and H given in [START_REF] Hu | Collaborative filtering for implicit feed-back datasets[END_REF] shall be applied several times to fully minimize the partial functions

W → D β (V|WH) and H → D β (V|WH)
. This leads to Algorithm 1, that we refer to as Block MM (BMM). Note that in common NMF practice, only one sub-iteration is used (L W = L H = 1), which still results in a descent algorithm thanks to the descent lemma (4).

Joint Majorization-Minimization

In contrast with the classic alternating scheme presented in Section 2.2, we develop in this section a joint MM (JMM) approach for solving Problem (1).

Construction of the auxiliary function

In order to apply the MM scheme, we start by looking for a suitable auxiliary given in Table 1.

Next, we majorize the convex and concave terms of d β separately. The methodology follows [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF], except that none of the two factors W or H is treated

β ⌣ d β (v f n | [WH] f n ) ⌢ d β (v f n | [WH] f n ) dβ (v f n ) ] -∞, 1[\{0} -v f n β-1 [WH] β-1 f n [WH] β f n β v β f n β(β-1) 0 v f n [WH] f n log [WH] f n -(log v f n + 1) [1, 2] d β (v f n | [WH] f n ) 0 0 ]2, +∞[ [WH] β f n β -v f n β-1 [WH] β-1 f n v β f n β(β-1)
as a fixed variable. The convex term

⌣ d β is majorized using Jensen's inequality ⌣ d β (v f n | [WH] f n ) = ⌣ d β v f n K k=1 w f k h kn = ⌣ d β v f n K k=1 λfnk w f k h kn λfnk ≤ K k=1 λfnk ⌣ d β v f n w f k h kn λfnk (7) def = ⌣ G f n (W, H | W, H) ,
where the coefficient λfnk is defined as λfnk = wfk hkn ṽfn and we denote Ṽ = W H with coefficients ṽfn . The concave term

⌢ d β is majorized using the tangent inequality ⌢ d β (v f n | [WH] f n ) ≤ ⌢ d β (v f n | ṽfn ) + ⌢ d ′ β (v f n | ṽfn ) [WH] f n -ṽfn def = ⌢ G f n (W, H| W, H) .
Finally the overall auxiliary function G is given by

G = F f =1 N n=1 ⌣ G f n + ⌢ G f n + dβ . (8)
By construction, it satisfies the tight majorization properties

G(W, H | W, H) ≥ D β (V | W H) G( W, H | W, H) = D β (V | W H) ,
which ensure the descent property in [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. The expression of G coincides with the joint auxiliary function derived in [23, Table 2

in Section 4 and Appendix

A] for a different purpose (see Section 3.3.4).

We stress out that the auxiliary function G is a tight joint majorization of D β . This is in contrast with the BMM approach of Section 2.2 where two separate auxiliary functions G W and G H are built for W → D β (V|WH) and

H → D β (V|WH)
, respectively. A central difference in our JMM approach is the definition of the coefficients { λfnk } which depend on both current iterates W and H and do not lead to a simplification of the term w f k h kn / λfnk in [START_REF] Fu | Nonnegative matrix factorization for signal and data analytics: Identifiability, algorithms, and applications[END_REF]. This is in contrast with BMM, where, say, W would be treated as a fixed variable and the term w f k h kn / λfnk simplifies to ṽfn h kn / hkn , allowing for closed-form minimization of G. Furthermore, the auxiliary function G is not jointly convex, due to the bilinear terms w f k h kn . It is however bi-convex, i.e., convex w.r.t W (resp., H) given H (resp., W).

Minimization step

The auxiliary function G does not appear to have a closed-form minimizer in W and H. Neither is it convex, which makes it difficult to minimize globally.

While minimizing G jointly in W and H is hard, we can still perform an alternating minimization on the matrices W and H which results in the following updates

W ←-W.    V ( W H) .(2-β) [χ 1,β (H, H)] ⊤ ( W H) .(β-1) [χ 2,β (H, H)] ⊤    .γ(β) H ←-H.    [χ 1,β (W, W)] ⊤ V ( W H) .(2-β) [χ 2,β (W, W)] ⊤ ( W H) .(β-1)    .γ(β) , (9)
where γ(β) is defined as in ( 6) and

χ 1,β (H, H) =    H.(2-β) H .(1-β) if β ≤ 2 H if β > 2 , χ 2,β (H, H) =    H if β < 1 H .β H.(β-1) if β ≥ 1 .
The updates [START_REF] Cichocki | Generalized Alpha-Beta divergences and their application to robust nonnegative matrix factorization[END_REF] are obtained by cancelling the partial gradients of G. Since G is not convex, the alternating minimization may only lead to a critical point instead of a global minimum. Nevertheless, in order to decrease the loss function D β , it is enough to decrease G thanks to the descent lemma (4). This is easily achieved by initializing the updates ( 9) with ( W, H), leading to Algorithm 2.

Note that, despite yielding a multiplicative update, the JMM updates given in ( 9) are conceptually different from the BMM ones: the JMM updates aim at minimizing G given W, H while the BMM updates aim at minimizing D β (V|WH) w.r.t H given W, or w.r.t W given H. Further comments and discussion are given in the next section.

Discussion

Special cases

For the values of notorious importance 0, 1 and 2 of β, the multiplicative updates (9) can be written in a simpler form as some of the exponents cancel out and make the corresponding terms vanish. Similar simplified updates are also available for the classic update rules [START_REF] Hu | Collaborative filtering for implicit feed-back datasets[END_REF]. These simplified formulae are shown in Table 2 for the factor matrix H, where the matrix 1 represents the matrix of dimension F × N whose all entries equal 1.

Computational advantages of JMM

The new update rules (9) have a similar structure to the classic multiplicative updates (5) with a notable difference regarding the matrices W and H. These matrices, named Wi and Hi in Algorithm 2, remain constant in the subiterations and allow for some computational savings w.r.t BMM when updating

β BMM JMM 0 H ← H. W ⊤ V (WH) .2 W ⊤ 1 WH . 1 2 H ← H. W.2 W ⊤ V ( W H) .2 W ⊤ 1 W H . 1 2 1 H ← H. W ⊤ V WH W ⊤ 1 H ← H. W⊤ V W H W ⊤ 1 2 H ← H. W ⊤ V W ⊤ WH H ← H. W ⊤ V W .2 W ⊤ W H
W l and H l . For instance, the computation of Ṽi = Wi Hi is performed only once per outer iteration in step 5 of Algorithm 2 whereas in Algorithm 1 the product WH has to be computed at each update of W l (product W l Ȟi at step 6) and at each update of W l (product Wi+1 H l at step 11). Our proposed update hence saves here a matrix product per sub-iteration. Similar computational savings can be obtained by storing for example, the entry-wise ratio of V and Ṽ in the update of JMM for β = 1.

Table 3 summarizes the computational savings (and extra divisions) of JMM w.r.t BMM for the different values of β at each iteration, i.e., for one update of Wi / Wi and one update of Hi / Ȟi . For a fair comparison, we pick L W and L H equal to L. The part multiplied by L in the expressions in Table 3 corresponds to the computational savings obtained at each sub-iteration for W l and H l while the remaining part corresponds to the extra cost of computing matrices that are constant through the L sub-iterations. We especially emphasize the values 0, 1, and 2 of β, which enjoy even larger computational savings thanks to the simplifications shown in Table 2. Since these three cases are the most common in practice, the update (9) brings a very welcome speedup compared with [START_REF] Hu | Collaborative filtering for implicit feed-back datasets[END_REF]. It turns out that the largest saving is in the case where β is equal to 0 or 1, i.e., the IS and the KL divergences. For other general values of β,

Multiplications

Divisions Additions

β = 0 (2L -1)F N K + 2LF N 2(L -1)F N -L(F K + KN ) (2L -1)F N (K -1) β = 1 (4L -3)F N K (2L -1)F N (4L -3)F N K -(L -1)(F N + F K + KN ) β = 2 (2L -1)F N K -L(F K + KN ) (2L -1)F N (K -1) β ∈]1, 2[ (2L -1)F N K -L(F K + KN ) -L(2F N -F K -KN ) (2L -1)F N (K -1) β > 2 (2L -1)(F N K + F N ) -L(F K + KN ) (2L -1)F N (K -1) β < 1 (2L -1)F N K -LK(F K + KN ) -F N (2L -1)F N (K -1)
the JMM updates incur extra divisions in ( 9), especially for β in ]1, 2[, which mitigate the computational savings.

Failure of heuristic updates

A common heuristic in BMM consists in setting γ(β) to 1 in the update [START_REF] Hu | Collaborative filtering for implicit feed-back datasets[END_REF], even for values of β outside [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF][START_REF] Smaragdis | Static and dynamic source separation using nonnegative factorizations: A unified view[END_REF]. Indeed, it has been empirically observed that this leads to an equally good factorization while decreasing the number of iterations to reach convergence. For values of β in [0, 1], the authors in [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF] have proved that this heuristic corresponds to a Majorization-Equalization scheme which produces larger steps than the MM method. Nevertheless, deriving theoretical support for this heuristic for other values of β is still an open problem.

Setting γ(β) = 1 for all β did not lead to similar findings for JMM. While we did observe that the objective function decreases at every iteration under the heuristic, worse solutions were obtained (i.e., corresponding to higher values of the objective function in general). This might be due to the non-convexity of the JMM auxiliary function G, in which case the Majorization-Equalization principle makes less sense.

Convergence of the iterates of JMM

In standard NMF practice, BMM is used with

L W = L H = 1 sub-iterations.
The convergence of the iterates of BMM in this setting can be proven for slightly modified NMF problems that essentially ensure the coercivity of the objective function on its domain (loosely speaking, f (θ) → ∞ whenever θ → ∞).

In [START_REF] Zhao | A unified convergence analysis of the multiplicative update algorithm for regularized nonnegative matrix factorization[END_REF] this is ensured by augmenting the objective function with an ℓ 1 regularization term on W and H. Then the convergence of the iterates can be invoked using the Block Successive Upper-bound Minimization (BSUM) framework of [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF]. In [START_REF] Takahashi | A unified global convergence analysis of multiplicative update rules for nonnegative matrix factorization[END_REF], coercivity is ensured by replacing the nonnegativity constraint with a strict positivity constraint W, H ≥ ε. Then the convergence of the iterates can be invoked using Zangwill's convergence theorem by astutely reformulating BMM (with L W = L H = 1) as follows: ∀l ≥ 1,

W l+1 = arg min W≥ε G(W, H l | W l , H l ) (10)
H l+1 = arg min H≥ε G(W l+1 , H | W l+1 , H l ). (11)
This is how the joint auxiliary function G given by ( 8) is introduced in [START_REF] Takahashi | A unified global convergence analysis of multiplicative update rules for nonnegative matrix factorization[END_REF],

namely as a convenient way to derive BMM from a unique function of four variables, rather than using separate auxiliary functions (of two variables) for each sub-problem. The fundamental difference between BMM and our novel JMM approach is that the second occurrence of W l+1 in ( 11) is left unchanged in JMM (i.e., it remains W l ). This seemingly trivial change has significant computational implications in practice (and can only be justified by the joint MM framework that we introduced).

The proofs of convergence in [START_REF] Zhao | A unified convergence analysis of the multiplicative update algorithm for regularized nonnegative matrix factorization[END_REF][START_REF] Takahashi | A unified global convergence analysis of multiplicative update rules for nonnegative matrix factorization[END_REF] heavily rely on the block-structure of BMM and the strict convexity of the auxiliary functions [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]

Experimental Results

We provide in this section extensive numerical comparisons of BMM and JMM using various datasets (face images, audio spectrograms, song play-counts, hyperspectral images) with diverse dimensions.

Set-up

Our implementation for JMM follows Algorithm 2 while the one for BMM follows Algorithm 1. Some additional practical considerations are detailed below.

Influence of the number of sub-iterations

Algorithm 2 dictates that we update W and H several times in the inner loop to fully minimize (W, H) → G(W, H| W, H) (without updating or recomputing the tilde matrices). Nevertheless, this does not turn out particularly advantageous in practice. Indeed, in all our simulations, we observe that only the first iteration of ( 9) results in a significant decrease of G and that the additional sub-iterations yield only negligible improvement. This is illustrated in Figure 1 which displays the objective function values obtained with JMM and BMM (as a function of the outer iterations) for different numbers of sub-iterations L, L W , and L H , using the Olivetti face dataset (see Section 4.2) and β = 1. Figure 1 shows that the plots for JMM with L = 1 and L = 10 sub-iterations are nearly overlapping. As such, like in traditional NMF practice, in the following we only use one sub-iteration of W and H in our implementations of JMM and BMM. This means we set L = 1 in Algorithm 2 and L W = L H = 1 in Algorithm 1. By doing so, the BMM and JMM updates of W coincide and only the update of H changes (with still a significant gain in performance). Note that we have observed empirically that inverting the order of the updates in [START_REF] Cichocki | Generalized Alpha-Beta divergences and their application to robust nonnegative matrix factorization[END_REF] does not have any impact on the number of iterations before convergence nor on the quality of the obtained solutions.

Initialization and stopping criterion

The initializations (W init , H init ) for BMM and JMM are drawn randomly according to a half-normal distribution. Our stopping criterion for both algorithms is based on the relative decrease in the objective function D β . More precisely, the algorithms are stopped when

D β (V | W H) -D β (V | WH) D β (V | WH) ≤ ǫ ,
where ǫ is a tolerance set to 10 -5 , W and H are the current outer-loop iterates while W and H are the previous ones (either W and Ȟ for BMM or W and H for JMM). Furthermore, in order to remove the scaling ambiguity inherent to NMF, we normalize W and H at the end of each outer-loop iteration for both BMM and JMM. The normalization consists in dividing each column of W by its ℓ 2 norm and scaling the rows of H accordingly.

Handling zero values and numerical stability

The β-divergence D β (x|y) may not be defined when either x or y takes zero values. This is for instance the case when β is set to 0 or 1 due to the quotient and the logarithm that appear. As such, we recommend minimizing

D β (V + κ1|WH + κ1
) with a small constant κ instead of D β (V|WH) for numerical stability. This first simply amounts to replacing V by V + κ1 in the previous derivations. Then, by treating κ as a (K + 1) th constant component like in [START_REF] Févotte | Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization[END_REF], we may simply replace the product W H (resp. WH) by W H + κ1 (resp. WH + κ1) in ( 9) (resp. ( 5)).

Simulation environment

All the simulations have been conducted in Matlab 2020a running on an Intel i7-8650U CPU with a clock cycle of 1.90GHz shipped with 16GB of memory. 1In each of the following experimental scenarii, we compare the factorization obtained by JMM and BMM from 25 different initializations (W init , H init ).

Performance evaluation

We compare BMM and JMM both in terms of computational efficiency (CPU time) and quality of the returned solutions. To assess the latter, we use the value of the normalized objective function

D β (V| W H) F N
at the solution ( W, H) returned by the NMF algorithms. We also consider the KKT residuals [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF] to measure the distance to the first order optimality conditions and attest whether the algorithms reached a critical point of the criterion D β . The KKT residuals are expressed by

res(W) = min{W, [(WH) .(β-2) .(WH -V)]H ⊤ } 1 F K res(H) = min{H, W ⊤ [(WH) .(β-2) .(WH -V)]} 1 KN .

Factorization of face images

In the context of image processing, NMF can be used to learn part-based features from a collection of images [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. The columns of the matrix V correspond to the vectorization of the different images. Besides, the factor W represents the dictionary of image features, and the matrix H contains the activation encodings.

We compare the BMM and JMM methods on a face images dataset, the Olivetti dataset from AT&T Laboratories Cambridge [START_REF] Samaria | Parameterisation of a stochastic model for human face identification[END_REF], which contains 400 greyscale images of faces with dimensions 64 × 64. The corresponding matrix V thereby has dimensions 4096 × 400. We set K to 10. We consider the values 0, 1, and 2 for the parameter β, which correspond respectively to IS divergence, KL divergence, and squared Euclidean distance (the latter two being the most common in image processing).

Figure 2 shows the computation times for β = 2 and β = 1, as well as the values of the normalized objective function produced by both BMM and JMM, for the 25 runs. The same random initialization is used by both methods. Note that we use a logarithmic scale for the CPU time and that the y-axis for the objective function does not start at zero. We observe that JMM is always faster than BMM while yielding solutions with a similar quality in terms of the objective D β . The corresponding average CPU time as well as the resulting acceleration ratios are given in Table 4. We notice that the computational saving is higher when using the IS divergence, which confirms our analysis from Section 3.3.2. Remark that we measure here the global CPU time and not the time per iteration: since the auxiliary function is different for BMM and JMM, the trajectory in the parameter space is also different and thus the two algorithms do not require the same number of iterations before convergence.

Nevertheless, we observe that the number of iterations for both algorithms has a similar order of magnitude. Since the iterations of JMM are cheaper, its global CPU time is lower. This remark also holds for higher dimensional dataset such as the one in Section 4.4.

Note that while MM algorithms are prevalent for β-NMF in general, many other algorithms have been designed for the specific case β = 2 (quadratic loss) [START_REF] Gillis | Nonnegative Matrix Factorization[END_REF]. In that case, some algorithms are notoriously more efficient than BMM, see, e.g., [START_REF] Gillis | Nonnegative Matrix Factorization[END_REF]Chapter 8.2]. Though JMM improves on BMM, it may not compete with these other algorithms. Still, JMM, like BMM, is free of hyperparameters and very easy to use, and remains a very convenient option for the general practitioner.

Finally, we computed the KKT residuals for W and H returned by JMM and BMM, for all runs and considered values of β. We observed that they range from 1.10 -5 to 1.10 -1 , indicating that both methods converge to a critical point of D β . As a matter of fact, an additional significant observation is that in all cases here, both JMM and BMM return the same solution ( Ŵ, Ĥ) up to permutation of their columns and up to some round-off errors. Nevertheless, the trajectory of the iterates may differ.

Factorization of a spectrogram

We now consider NMF of audio magnitude spectrograms. In this context the factor W contains elementary audio spectra with temporal activations given by H [START_REF] Smaragdis | Static and dynamic source separation using nonnegative factorizations: A unified view[END_REF].

We generate the spectrogram of an excerpt from the original recording of the song "Four on Six" by Wes Montgomery. The signal is 50-seconds long with a sampling rate of 44.1 kHz. The spectrogram is computed with a Hamming window of length 2048 (46ms) and an overlap of 50%. This results in a data matrix V of dimensions 1025 × 2152.

In this section, we set β = 0, which is a common value in audio signal processing [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF], and set K to 10. The results obtained with BMM and JMM are shown on Figure 3. The average computation time for BMM is 291s while the one for JMM is 41s, which yields an average acceleration of 86%. We observe a very substantial benefit of JMM over BMM in terms of CPU time in this case.

This confirms our conclusion from Section 3.3.2 that JMM reaches its highest potential with NMF based on the IS divergence. We observe here again that JMM and BMM return the same solution ( Ŵ, Ĥ) up to permutation of their columns.

Factorization of song play-counts

NMF may be used in recommendation systems based on implicit user data.

In this case, the matrix V contains information about the interactions of users with a collection of items. The factor W may extract user preferences while H represents item attributes, see, e.g., [START_REF] Hu | Collaborative filtering for implicit feed-back datasets[END_REF].

We here consider the TasteProfile dataset [START_REF] Bertin-Mahieux | The million song dataset[END_REF] which contains counts of songs played by users (e.g., of a music streaming service). Similarly to [START_REF] Gouvert | Ordinal non-negative matrix factorization for recommendation[END_REF] and many other papers using this dataset, we apply a preprocessing to the original data to keep only users and songs with more than a given number of interactions (here set to 20). This results in a large and sparse data matrix V of dimensions 16 301 × 12 118 with about 0.6% nonzero values. We set β to 1 (a common choice in recommender systems, because its corresponds to the log-likelihood of a Poisson model that is natural for count data), and K = 50.

Comparative results are displayed in Figure 4. We observe that on average, 20 Similarly to our observations of Section 4.2, we observe that, while BMM has generally a shorter trajectory than BMM, the number of iterations of both methods has the same order of magnitude. Since the cost of BMM per iteration is higher than JMM, this results in a higher overall CPU time.

Factorization of hyperspectral images

A hyperspectral image is a multi-band image that can be represented by a nonnegative matrix: each row represents a spectral band while each column represents a pixel of the image. Applying NMF to such matrix data allows extracting a collection of individual spectra representing the different materials arranged in the matrix W, as well as their relative proportions given by the matrix H, see, e.g., [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF].

We consider a hyperspectral image acquired over Moffett Field in 1997 by the Airborne Visible Infrared Imaging Spectrometer [START_REF]Aviris free data[END_REF]. The image contains 50×50 pixels over 189 spectral bands, which leads to a matrix V of dimensions 189 × 2500. We consider β = 2 and β = 1.5. The latter value in particular was shown to be an interesting trade-off between Poisson (β = 1) and additive Gaussian (β = 2) assumptions for predicting missing values in incomplete versions of this dataset, see [START_REF] Févotte | Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization[END_REF]. The factorization rank K is set to 3, a standard choice with this dataset (extraction of vegetation, soil and water).

Comparative results are given in Figure 5. On the top figure corresponding to β = 2, we observe that JMM is faster than BMM with an acceleration of 31% on average. On the other hand, for β = 1.5, we observe that BMM is faster in general (though not always). This confirms again our analysis in Section 3.3.2:

when β is different from 0, 1, and 2, the benefit of JMM may be counterbalanced by the additional divisions in the general formulae (9).

Summary of the experiments

Table 4 summarizes the average and the standard deviation of the CPU times for all the datasets considered. It can be observed that JMM yields a significant speed-up for the widely used settings when β equals to 0, 1, or 2.

This acceleration is observed on small, medium and large size datasets. The corresponding 95% confidence intervals are given in Table 5.

Conclusion

In this paper, we have presented a joint MM method for NMF with the turned by the classic block MM scheme. As a matter of fact, we have noted that both algorithms usually return the same solutions up to column permutation.

The computational efficiency of the proposed updates has been demonstrated on datasets with diverse characteristics. In future work, we intend to study the convergence of the iterates of JMM. This is a challenging topic because of the non-convexity of the majorizer which underpins our approach. Such a difficulty explains the scarcity of results in the literature on the convergence of MM algorithms with non-convex auxiliary functions. Another promising topic would in designing stochastic versions of JMM for the factorization of massive datasets [START_REF] Mensch | Stochastic subsampling for factorizing huge matrices[END_REF][START_REF] Pu | Stochastic mirror descent for lowrank tensor decomposition under non-Euclidean losses[END_REF]. Initialize W 1 to Wi 5:

for l = 1 . . . L W do 6:

Update W using ( 5)

W l+1 ←-W l .   (W l Ȟi ) .(β-2) .V Ȟ⊤ i (W l Ȟi ) .(β-1) Ȟ⊤ i   .γ(β) 7:
end for 8:

Wi+1 ← W LW +1

9:

Initialize H 1 to Ȟi Ṽi ← Wi Hi

5:

Initialize (W 1 , H 1 ) to ( Wi , Hi )

6:

for l = 1 . . . L do 7:

Update W using ( 9)

W l+1 ← Wi .    V Ṽ.(2-β)

  function G : R F ×K + ×R K×N + → R + for the function (W, H) → D β (V|WH). We observe in (2) that D β is a sum of F N β-divergences between scalars. Our approach is to construct an auxiliary function for each summand. Following [13], the β-divergence d β , taken as a function of its second argument, can be decomposed into the sum of a convex term ⌣ d β , a concave term ⌢ d β , and a constant term dβ . The definitions of these three terms for the different values of β are

Figure 1 :

 1 Figure 1: Impact of the number of sub-iterations in the minimization step of Algorithm 2.The box is a zoom-in on iterations 750 to 881. In the latter, BMM (L W = L H = 10) and JMM (L = 10) are overlapping.

  Using KL divergence (β = 1).

  Using quadratic loss (β = 2).

Figure 2 :

 2 Figure 2: Comparative performance using the Olivetti dataset (K = 10).

Figure 3 :

 3 Figure 3: Comparative performance with a spectrogram (K = 10, β = 0).

Figure 4 :

 4 Figure 4: Comparative performance with TasteProfile (K = 50, β = 1).

β

  -divergence. Our algorithm relies on the alternating minimization of a nonconvex auxiliary function and leads to new multiplicative updates of the factor matrices. These new updates are variants of the classic multiplicative updates and are equally simple to implement. They can lead to a significant computational speedup, especially for the Itakura-Saito and Kullback Leibler divergences, and the quadratic loss. Fortunately, the three later cases are the three 22

Figure 5 :

 5 Figure 5: Comparative performance with Moffet hyperspectral image (K = 3).

Algorithm 1 BMM 1 :i to 1 2:

 111 Input: Nonnegative matrix V and initialization (W init , H init ) Output: Nonnegative matrices W and H such that V ≈ WH Initialize Initialize ( Wi , Ȟi ) to (W init , H init )

i [χ 1 ( 13 :

 113 ,β (H l , Hi )] ⊤ Ṽ.(β-1) i [χ 2,β (H l , Hi )] ,β (W l+1 , Wi )] ⊤ V Ṽ.(2-β) i [χ 2,β (W l+1 , Wi )] ⊤ Ṽ.(β-1) Wi+1 , Hi+1 ) ← (W L+1 , H L+1 ) return ( Wi , Hi )

Table 1 :

 1 Decomposition of d β for the different values of β.

Table 2 :

 2 Simplified updates of H for β = 0, 1, 2.

Table 3 :

 3 Computational savings and extra divisions brought by JMM compared to BMM. The table reports the difference between the number of operations required by JMM and BMM per iteration of the outer loop, using L W = L H = L. Benefits of JMM are highlighted in bold font.

Table 4 :

 4 Summary for the statistics on CPU times and the acceleration in all the tested dataset. The displayed values are the averages while the standard deviations are indicated within parenthesis.

			BMM	JMM	Acc.
	Olivetti	(β = 2)	55s (8s)	36s (6s)	35%
	Olivetti	(β = 1)	40s (9s)	34s (6s)	16%
	Olivetti	(β = 0)	229s (32s)	63s (8s)	72%
	Spectrogram (β = 0)	291s (90s)	41s (8s)	86%
	TasteProfile (β = 1) 6776s (439s) 5909s (492s) 13%
	Moffet	(β = 2)	39s (13s)	24s (5s)	35%
	Moffet	(β = 1.5)	107s (22s)	123s (26s)	-18%

  Algorithm 2 JMM Input: Nonnegative matrix V and initialization (W init , H init ) Output: Nonnegative matrices W and H such that V ≈ WH 1: Initialize i to 1 2: Initialize ( Wi , Hi ) to (W init , H init )

	11:	Update H using (5)
		H l+1 ←-H l .	 	W⊤ i+1 ( Wi+1 H l ) W⊤ i+1 ( Wi+1 H l ) .(β-1) .(β-2) .V	 	.γ(β)
	12:	end for
	13:	Ȟi+1 ← H LH +1
	14:	Increment i
	15: until Convergence

10: for l = 1 . . . L H do 16: return ( Wi , Ȟi ) 3: repeat 4:

Matlab code is available at https://arthurmarmin.github.io/research.html.

Acknowledgment

The authors acknowledge Rémi Flamary, Jérôme Idier, Paul Magron and Emmanuel Soubies for discussions related to this work.

This work is supported by the European Research Council (ERC FACTORY-CoG-6681839), the French Agence Nationale de la Recherche (ANITI, ANR-19-P3IA-0004) and the National Research Foundation, Prime Minister's Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) programme.


Kien Lam Ung

Chloé Poussineau

Julie Couston

B Alsarraf

Mickaël Blaise
email: mickael.blaise@irim.cnrs.fr

Husam M A B Alsarraf

Crystal structure of MAB_4123 a putative flavin-dependent monooxygenase from Mycobacterium abscessus

Keywords: Mycobacterium abscessus, organosulfur compounds, two-component flavindependent monooxygenase

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1.Introduction

Monooxygenase proteins use flavins as cofactors to catalyze the breakage of the bond between two oxygen atoms. This reaction results in the addition of one oxygen atom into the substrate and the reduction of the second oxygen into water (Ellis, 2010). In this large enzyme family, the two-component flavin-dependent enzymes are only described in bacteria. The twocomponent flavin-dependent monooxygenases with the help of either flavin adenine dinucleotide (FAD) or flavin mononucleotide (FMN) oxidize compounds found in the bacterial environment.

These enzymes capable of desulfurization of sulfonated compounds (Fig. 1) are of particular interest to biotechnology as they can notably degrade sulfur-compounds from fossil fuels, a source of environmental pollution.

DszC from Rhodococcus sp. XP (Liu et al., 2014) or TdsC from Paenibacillus sp. A11-2 (Hino et al., 2017) catalyze the two first steps of dibenzothiophene (DBT) degradation (Fig. 1) leading to the formation of DBT sulfone that is further processed by the dibenzothiophene sulfone monooxygenase DszA/TdsA and the 2-(2'-hydroxyphenyl)-benzene sulfinate (HPBSi) desulfinase DszB/TdsB (Piddington et al., 1995;Ishii et al., 2000;Kirimura et al., 2004;Ohshiro et al., 2005;Konishi et al., 1997). In this pathway, the DszC and TdsC proteins have been well characterized biochemically and structurally, hence enabling the understanding of their catalytic function at a molecular level (Liu et al., 2014;Hino et al., 2017).

Similarly, the two-component flavin-dependent monooxygenase MsuC from Pseudomonas fluorescens was shown to catalyze the second step of the dimethylsulfone (DMSO 2 ) degradation pathway (Kertesz et al., 1999;Wicht, 2016). MsuC converts methanesulfinate (MSI -) to methanesulfonate (MS -) (Fig. 1) which is finally processed to sulfite by MsuD that can be assimilated by different bacterial synthesis pathways (Soule et al., 2020). Crystallographic investigation of MsuC confirmed the close structural similarity with the DszC/TszC enzymes (Soule et al., 2020).

As exemplified above, numerous bacteria possess these organosulfur detoxification pathways enabling organosulfur compounds degradation.

Pathogenic bacteria such as

Mycobacterium goodii a rapidly growing nontuberculous mycobacterium and an opportunistic human pathogen can possess as well this pathway. Of interest, the DszC enzyme from Mycobacterium goodii X 7 B (Li et al., 2009) processes DBT using either FMNH 2 or FADH 2 as a cofactor.

So far, the mycobacterial two-component flavin-dependent monooxygenases have not been investigated from a structural point of view. To fill this gap, we investigated the putative twocomponent flavin-dependent monooxygenase encoded by MAB_4123 from Mycobacterium abscessus, an emerging human pathogen, highly resistant to antibiotherapy (Johansen et al., 2020).

We report herein the crystal structure of MAB_4123 at high resolution.

2.Materials and Methods

Cloning, expression and purification of MAB_4123

The wild-type DNA sequence of MAB_4123 (Uniprot: B1MIH5) was amplified from the genomic DNA of M. abscessus reference strain 19977 using a pair of primers described in Table I.

The resulting PCR fragment was cloned into pET32(a) expression plasmid between the KpnI-XhoI restriction sites. The resulting plasmid encoded for a fusion protein with a thioredoxin tag, hexahistidine, thrombin site, S-tag, and Tobacco Etch Virus (TEV) protease cleavage site at the Nterminus of MAB_4123. This plasmid was transformed into E. coli BL21 strain (New England Biolabs, C2527) strain containing the pRARE2 plasmid and used for protein production.

Transformed colonies were inoculated into LB medium supplemented with 200 μg.mL -1 ampicillin and 30 μg.mL -1 chloramphenicol and grown overnight at 37 °C under agitation.

Crystallization

Prior to crystallization, 2 mM dibenzothiophene and 4% dimethylsulfoxide (DMSO) were added to the protein solution concentrated at 9.5 mg.mL -1 in 50 mM Tris pH 8.0 and 0.2 M NaCl.

Crystallization drops were set-up in 96 well MRC plates mixing 0.8 μl of protein solution with 0.8 μl of the reservoir solution. Crystals were obtained at 291 K in 1.5 M Ammonium phosphate dibasic and 0.1 M Tris pH 8.5 and cryo-protected in 1.5 M Ammonium phosphate dibasic, 0.1 M Tris pH 8.5 and 20% glycerol (v/v) before being cryo-cooled in liquid nitrogen.

Table 2 Crystallization

Method

Data collection and structure solving

The X-ray dataset was collected on the ID30A1-MASSIF1 beamline (Svensson et al., 2015) at the European Synchrotron Radiation Facility, Grenoble, France. The dataset was recorded on a PILATUS3-2M detector (Dectris) at a wavelength of 0.965 Å (12.842 keV) and a crystal-to-detector distance of 230.6 mm. A total of 3160 frames were collected with an exposure time of 0.11 s, a rotation range of 0.05° and full beam transmission. Data were processed, scaled, and merged with XDS (Kabsch, 2010). The structure was solved by molecular replacement using the Molrep (Vagin & Teplyakov, 2010) program with the structure of the dibenzothiophene monooxygenase (DszC)

from Rhodococcus erythropolis (PDB id 4JEK) as a search model and sharing 37% protein sequence identity with MAB_4123. The tetramer obtained from the first molecular replacement was used as a search model to perform an additional molecular replacement step using Phaser (McCoy et al., 2007). The model was then automatically adjusted with phenix.autobuild from the Phenix software suite (Liebschner et al., 2019) followed by several cycles of manual adjustment and refinement with Coot (Emsley & Cowtan, 2004) and phenix.refine (Liebschner et al., 2019) respectively. The statistics for the data collection and structure refinement are reported in Tables 3 and4.

Ligand in silico docking

FMN was first subjected to molecular dynamics energy minimization by Merck molecular force field (MMFF94) using the PyRx software suite (Dallakyan & Olson, 2015). Hydrogen atoms were added both to protein and ligand. The docking was performed with AutoDock Vina (Trott & Olson, 2010) run in PyRx using the structure of the MAB_4123 dimer as a receptor and with residue

Arg361 defined as flexible. The search was performed with a docking box centered as followed :

x=-65.8, y=-92.4, and z=-73.8 with dimension (Å) X=25.9, Y=30.1 and Z=39.7, and exhaustiveness defined to 8.

3.Results and Discussion

Purification, crystallization and structure determination of MAB_4123

The full-length MAB_4123 gene was cloned into pET-32a and expressed as a fusion protein with the thioredoxin and multi-histidine tags in E. coli. The protein was purified through three steps of purification involving two steps on IMAC column and one purification step on size-exclusion chromatography (SEC) (Fig. 2). We could crystallize the protein in presence of a potential ligand, dibenzothiophene (DBT) diluted in dimethyl sulfoxide (DMSO). These crystals belonging to the monoclinic P2 1 space group diffracted to a resolution of 2.1 Å. We solved the phase problem by molecular replacement using the structure of the DszC protein from Rhodococcus erythropolis (PDB id 4JEK). The asymmetric unit contains 16 monomers (Fig. 3a). We could model most of the residues apart from the first 11 residues in the N-terminus (Table I) and the last seventeen residues in the C-terminus as well as residues 136-139 and 174-175 in chains H and N respectively. In addition, we placed several phosphate ions found at the interface of two monomers. Packing analysis suggested the existence of four stable tetramers within the asymmetric unit (Figs. 3a and3b). The homotetramer assembly was further confirmed in solution by SEC (Fig. 2a). Indeed, the MAB_4123 elution peak on SEC corresponds to an apparent molecular weight of 162 kDa while the theoretical mass of the protein is 46 kDa. Such assembly was also attested in MAB_4123 close homologs notably TsdC (Hino et al., 2017) and MsuC (Soule et al., 2020).

MAB_4123 structural analysis

The protein adopts a typical acyl-CoA dehydrogenase fold (Fig. 3). The N-terminal αdomain (residues 10-129) is composed of five α-helices followed by the central domain (residues 130-232) made of two antiparallel β-sheets forming a β-sandwich and the C-terminal α-domain (residues 233-405) comprising five α-helices.

Using the ENDscript server (Robert & Gouet, 2014), we retrieved close homologous structures to MAB_4123. As stated in the introduction, DszC and TdsC were shown to be dibenzothiophene monooxygenases while MusC was recently demonstrated to be involved in the conversion of methanesulfinate (MSI -) to methanesulfonate (MS -) (Fig. 1). These three enzymes present respectively 35%, 34%, and 31% protein sequence identity with MAB_4123 (Fig. 4a). The Pseudomonas fluorescens (PDB id: 6UUG) (Fig. 4b). Overall, the structures are very similar with only a few noticeable differences. Firstly, the region between α6 and β7 ranging from residue 217 to 233, folded as a loop and two 3 10 helices in MAB_4123, is very similar to the one of MsuC but differs from the ones of DsbC/TdsC which is more folded with the presence of one α-helix and one 3 10 helix. The second noticeable difference lies in the length of loop 1 spanning from residues 18CDA31 to 137 which is of similar length in MsuC but slightly longer in DsbC/TdsC as twelve residues composed loop 1 in these proteins. It is also worth mentioning that the electron density map for loop 1 between strands β1 and β2 and loop 2 (271-286) between helices α6 and α7 is well defined in most MAB_4123 monomers (except for loop 1 in chain H) and in a closed conformation, these characteristics are more related to ligands/cofactor bound structures as proposed for the TdsC protein (Hino et al., 2017). The role of loop 2 is not well defined in these enzymes but loop 1 is near the active site residues. A third difference is the length of the loop separating strands β4 and β5 made of five residues in MAB_4123 while in TdsC/DsbC this loop has 4 additional residues.

Structures superimposition attests that MAB_4123 is a very close homolog of the DsbC/TdsC/MsuC enzymes, thus strongly suggesting that MAB_4123 belongs also to the flavindependent monooxygenase class of enzymes.

Cofactor and putative ligand binding sites.

To further investigate the potential function of MAB_4123, we assessed if the protein as its close homologs, could bind flavin cofactor and notably FMN. We tried to obtain the structure of MAB_4123 bound to FMN either by co-crystallization or soaking experiment, but without any success. Alternatively, we performed in silico docking of FMN into MAB_4123 with a search grid that was centered on one monomer and extended partly on the second monomer (Fig. 5a). The search volume encompassed the putative binding site that was defined by structural comparison with the TdsC structure bound to FMN (Hino et al., 2017). This approach led to identifying a potential binding site of FMN in MAB_4123 at the interface of two monomers as seen in its close homologs. In the best docking pose (ΔGº=-10.3 kcal.mol -1 ) sixteen residues form a binding cavity for FMN (Fig. 5b). Most of these amino acids are strictly conserved and semi-conserved in all MAB_4123 homologs (Figs. 3b and5d). Residues from one monomer mainly interact with the flavin moiety of FMN while the recognition of the phosphate group is mediated by two Arginines from two different monomers namely Arg137 and Arg361. In agreement with this docking pose, the phosphate group of FMN contacted by Arg137 and Arg361 from each monomer is positioned very similarly to a phosphate ion we could unambiguously identify and model (Fig. 5c). Indeed, the extra electron density in the vicinity of these two Arg is well-defined.

Additionally, we noticed two positive electron density peaks in a cavity formed by residues Ser102, Ser245, Trp106, Tyr99, His379, Ser134, and Ser132 (Fig. 5d). Interestingly, this unmodelled electron density is situated above the flavin ring of FMN and matches the position of DBT as seen in the cocrystal structure of TdsC bound to DBT and FMN (PDB id 5XDE) (Hino et al., 2017). There is good conservation of the position of the substrates and residues surrounding these extra electron density peaks between the MAB_4123 and TdsC structures so these two additional electron density peaks could be potential substrates of MAB_4123. However, we rule out that this could be DBT (which was present in the crystallization drop) as the electron density is discontinued and corresponds more to two small molecules. On top of DBT presence in the crystallization drop, 4% DMSO was also present in the crystallization condition. As it was shown that MsuC, a close homolog of MAB_4123, converts MSI -to MS -. MSI -is chemically related to DMSO and it is tempting to propose that this density could correspond to DMSO, at least for the extra density near Ser245 and Ser102 (Fig. 5d). However, as we could not obtain an exploitable dataset of crystal grown without DMSO, as we cannot guarantee that these positive electron density peaks are DMSO, we, therefore, left them unmodelled. Nonetheless, these observations suggest that MAB_4123 like MsuC could detoxify methanesulfinate or closely related chemical analogs.

Concluding remarks

To conclude, this study demonstrates the structural proximity of MAB_4123 to the TdsC/DszC and MsuC enzymes that are all members of the two-component flavin-dependent monooxygenase family and involved in the detoxification of organosulfur compounds. The putative role of MAB_4123 in a such pathway is also very well supported by the genomic organization of the M. abscessus genome. Indeed, MAB_4123 is part of a predicted operon where a potential nitrilotriacetate monooxygenase is encoded by MAB_4125. This protein harbors a bacterial luciferase predicted domain seen in the dibenzothiophene desulfurization enzymes A that is working in pair with DsbC in the dibenzothiophene degradation pathway. Further, several papers have reported dibenzothiophene degradation pathway in mycobacteria (Nomura et al., 2005) and notably the study on Mycobacterium goodii X7B describing the DsbC enzyme, which shares 37% sequence identity with MAB_4123 (Li et al., 2009). Further, blast search with the MAB_4123 sequence against the M. abscessus genome retrieved only the MAB_1422c protein as a paralog which displays only 29% identity with MAB_4123 and moderate sequence identity (<25%) with the DsbC/TdsC and MsuC enzymes. Altogether, it is highly probable, like other nontuberculous mycobacteria that M. abscessus possesses an organosulfur detoxification pathway where MAB_4123 is part of it. However, confirmation of MAB_4123 function awaits biochemical assay.

  Vapor drop diffusion method in sitting drops Plate type Swissci 96-Well 2-Drop MRC Crystallization Plates Temperature (K) 291 Protein concentration 9.5 mg.mL -1 Buffer composition of protein solution 50 mM Tris-HCl pH 8.0, 200 mM NaCl Composition of reservoir solution 0.1 M Tris-HCl pH 8.5 and 1.5 M Ammonium phosphate dibasic Volume and ratio of drop 0.8 μL+ 0.8 μL Volume of reservoir 50 μL

  MAB_4123 monomer displays an r.m.s.d. of 1.43 Å over 355 Cα with DszC from Rhodococcus erythropolis (PDB id 4JEK), an r.m.s.d. of 1.37 Å over 350 Cα with TdsC from Paenibacillus sp.A11-2 (PDB id: 5XB8) and finally an r.m.s.d. of 1.45 Å over 358 Cα with MsuC from

Figure 2

 2 Figure 2 Purification and oligomeric state of MAB_4123 a) Size-exclusion chromatogram from the Superdex 200 Increase 10/300 GL column. The void volume (8.58 mL) was determined with dextran blue (black arrow) and the calibration curve (black)was determined using β-amylase (Ve=11.98 mL, 200 kDa), bovine serum albumin (Ve=16.58 mL, 66 kDa), carbonic anhydrase (Ve=14.03 mL, 29 kDa), and cytochrome C (Ve=14.03 mL, 18.14 kDa). Proteins standard (250 to 500 μg injected two by two) and MAB_4123 were eluted in 20 mM Tris pH 8, 0.15 M NaCl and 10% glycerol at a flow-rate of 0.35 mL.min -1 . Kav was calculated as follows: Ve-Vo/Vt-Vo, where Ve is the elution of the protein, Vo is the void volume and Vt is the total volume of the column. The calibration curve was obtained by plotting the Kav against the logarithm (Log) of the molecular weight (M.W.) of each protein standard. MAB_4123 (red curve) eluted at 12.33 mL enabled the estimation of an apparent molecular weight in solution of 162 kDa. b) Coomassie-stained SDS-PAGE attesting of the high protein purity after the three steps purification procedure. 10 μg of protein were loaded on the gel.

Figure 3

 3 Figure 3 Overall crystal structure of MAB_4123 a) Depiction of the asymmetric unit composed of 16 monomers of MAB_4123 arranged into 4 tetramers. b) Overall structure of the MAB_4123 homotetramer shown as a cartoon representation and in two different views. The interface of interaction between two monomers is mainly formed by the C-terminal α-domain. The spheres indicate phosphate ions. c) Cartoon representation of the MAB_4123 monomer in two different views. The N-terminal α-domain is shown in blue, the central domain containing the antiparallel β-sheet in raspberry and the C-terminal α-domain in black. Loops 1 and 2 are indicated and are well-ordered.

Figure 4

 4 Figure 4 Structural comparison of MAB_4123 with close homologs a) Multiple sequence alignment based on structural homologs of MAB_4123 and performed with the ENDscript server (https://endscript.ibcp.fr) (Robert & Gouet, 2014). Only the three closest structures to MAB_4123 namely DszC from Rhodococcus erythropolis (PDB id 4JEK), TdsC from Paenibacillus sp.A11-2 (PDB id: 5XB8) and MsuC from Pseudomonas fluorescens (PDB id: 6UUG) are displayed. Strictly conserved and semi-conserved residues are coloured in red and yellow respectively. Secondary structures of MAB_4123 appear on top of the alignment. α, designs

Figure 5

 5 Figure 5 Identification of MAB_4123 cofactor and substrates binding sites a) Representation of the model used to performed in silico docking experiment. The search volume involved the MAB_4123 dimer and used to dock FMN is represented by the grid box in dark. The putative active site location is indicated with the red circle.

Table I Macromolecule production information

	Source organism Mycobacterium abscessus ATCC 19977
	DNA source: genomic DNA from Mycobacterium abscessus ATCC 19977
	Forward primer 5'-CTCACGGTACCGAGAACCTGTACTTCCAGGGTCATatgaccaccgccaccacctc-3' (letters in
	bold indicate the KpnI restriction site and in italic the sequence encoding the Tobacco Etch Virus protease
	cleavage site)
	Reverse primer 5'-GGATGCTCGAGtcacggcttcttgccctgtg-3' (letters in bold indicate the XhoI restriction
	site)
	Expression vector pET-32(a)
	Expression host E. coli BL21 (DE3) transformed with the pRARE2 plasmid
	The complete amino acid sequence of the construct after TEV cleavage (in bold extra amino acids from
	cloning strategy and TEV cleavage)
	GHMTTATTSPVHAYQGVSDTEFSEWEQVAARVAGELSATALTRDRANQNPIAEIELLRRYGLLSFAT
	AREFGGAGGSLVQALQLGRIIAAADGSIGQLLVYHYSNGVWTYILGSPTQREYISRGVGGHGWFQ
	GSVSNPRDPGITVTRTEEGYRVNGKRTFATGVAVADLITVLLYEAEPINAIIPSERDGLRFNDDWDNL
	GQRLTASGSVEFDNVLLRHDEVLTGLDEYSGLDGSRERRDGLRALFSQLIFVHLYLGIAEGALAAG
	VAYIRDKGRPWPEAHSTDVTEDPYHQQLLGRLSAGIAAGVALADSATKEFEQALAFGEAPTEAQW
	GALAIRVDQAKSVATEISLDVTHNIYQATGARSTANSVGLDIYWRNARTHTTHDPLPYRQREIGRH
	LLTDQWPSPRSLNGLGRLAETTQGKKP

Table 3 Data collection and processing

 3

	Diffraction source	ESRF, ID30A-1
	Wavelength (Å)	0.965
	Temperature (K)	100
	Detector	PILATUS3-2M
	Crystal-detector distance (mm)	230.6
	Rotation range per image (°)	0.05
	Total rotation range (°)	158
	Space group	P2 1
	a, b, c (Å)	131.91, 134.18, 206.80
	α, β, γ (°)	90, 90.68, 90
	Mosaicity (°)	0.082
	Resolution range (Å)	100-2.1 (2.2-2.1)
	Total No. of reflections	1238477 (114732)
	No. of unique reflections	413629 (41003)
	Completeness (%)	98.75 (98.3)
	Multiplicity	3 (2.8)

Values in parentheses are for the outer shell.

Table 4 Structure refinement

 4

	Values in parentheses are for the outer shell.
	Resolution range (Å)	47.94 -2.1 (2.175 -2.1)
	Reflections used in refinement	413476 (40988)
	Reflections used for R free	2016 (199)
	R work	0.178 (0.314)
	R free	0.212 (0.336)
	Number of non-hydrogen atoms	51940
	macromolecules	48833
	solvent	2986
	ligand	121
	Protein residues	6312
	R.m.s.d., bond lengths (Å)	0.002
	R.m.s.d., angles (°)	0.43
	Ramachandran favored (%)	98.60
	Ramachandran allowed (%)	1.40
	Ramachandran outliers (%)	0
	Rotamer outliers (%)	0.69
	Clashscore	2.76
	Average B-factor (Å 2 )	50.40
	Macromolecule	50.54
	Solvent	47.76

Acknowledgments

We thank the staff at the MASSIF-1 beamline at ESRF for support during data collection as well as the staff at the PXI and PXIII beamlines at the Swiss Light Source. We thank Dr. L. Kremer and Dr.

K. Brodolin for their support at the early stage and throughout the project, respectively. KLU PhD fellowship was supported by the National Research Agency (ANR-17-CE11-0008-01 -MyTraM to MB) and the Lundbeck Foundation, Denmark (R303-2018-2964 to HA). JC's PhD fellowship is financed by CNRS