Dawen Wu
email: dawen.wu@centralesupelec.fr

Abdel Lisser
email: abdel.lisser@l2s.centralesupelec.fr

Solving Constrained Pseudoconvex Optimization Problems with Deep Learning-Based Neurodynamic Optimization

Keywords: Constrained Pseudoconvex Optimization Problems, Neurodynamic Optimization, Neural Networks, Ordinary Differential Equations

In this paper, we consider Constrained Pseudoconvex Nonsmooth Optimization Problems (CPNOPs), which are a class of nonconvex optimization problems. Due to their nonconvexity, classical convex optimization algorithms are unable to solve them, while existing methods, i.e., numerical integration methods, are inadequate in terms of computational performance. In this paper, we propose a novel approach for solving CPNOPs that combines neurodynamic optimization with deep learning. We construct an initial value problem (IVP) involving a system of ordinary differential equations for a CPNOP and use a surrogate model based on a neural network to approximate the IVP. Our approach transforms the CPNOP into a neural network training problem, leveraging the power of deep learning infrastructure to improve computational performance and eliminate the need for traditional optimization solvers. Our experimental results show that our approach is superior to numerical integration methods in terms of both solution quality and computational efficiency.

Introduction

Constrained nonlinear optimization problems involve finding the best solution among a set of possible solutions by minimizing or maximizing an objective function. These problems are prevalent in various fields such as engineering, physics, finance, and management, with a wide range of applications. They can be divided into two groups based on the nature of the objective or constraint functions: convex and nonconvex optimization problems. Convex optimization problems, which include linear programming and quadratic programming, are a special class of nonconvex optimization problems and have been studied extensively.

Methods such as the primal-dual interior point method have been developed to solve them efficiently [START_REF] Boyd | Convex optimization[END_REF][START_REF] Nocedal | Numerical optimization[END_REF]. Nonconvex optimization problems, however, are more complex and commonly solved through gradient descent-based algorithms, which often struggle to converge to the global optimal solution [START_REF] Jain | Non-convex optimization for machine learning[END_REF][START_REF] Jiang | Structured nonconvex and nonsmooth optimization: algorithms and iteration complexity analysis[END_REF].

In this paper, we focus on a specific type of constrained nonconvex optimization problem known as constrained pseudoconvex nonsmooth optimization problems (CPNOPs). A CPNOP has an objective function that is both pseudoconvex, meaning that it is not strictly convex, and nonsmooth, meaning that it is not everywhere differentiable. Due to the pseudoconvexity and non-smoothness of the objective function, traditional convex optimization algorithms are not applicable to solve it.

CPNOPs are typically solved by neurodynamic optimization, which involves constructing an initial value problem (IVP) consisting of a system of ordinary differential equations (ODEs) [START_REF] Liu | A one-layer recurrent neural network for constrained pseudoconvex optimization and its application for dynamic portfolio optimization[END_REF][START_REF] Xu | A neurodynamic approach to nonsmooth constrained pseudoconvex optimization problem[END_REF][START_REF] Liu | A one-layer recurrent neural network for nonsmooth pseudoconvex optimization with quasiconvex inequality and affine equality constraints[END_REF]. The final state in the solution of this IVP represents the predicted solution of the CPNOP. These ODE systems are highly nonlinear and have no analytical solutions. Therefore, numerical integration methods, such as the Runge-Kutta method [START_REF] Dormand | A family of embedded runge-kutta formulae[END_REF] or the backward differentiation formula [START_REF] Shampine | The matlab ode suite[END_REF], are often used to approximate their solutions numerically. Motivation. However, numerical integration methods can be computationally expensive and can produce inaccurate solutions, particularly if the ODE system derived from the CPNOP is stiff. In addition, these methods are not well suited to solving CPNOPs as they require the calculation of all previous states to reach the desired final state, making them inefficient. Therefore, a more efficient and accurate approach is needed to solve CPNOPs.

Related Works

Neuraodynamic Optimization

Neurodynamic optimization is a class of methods that model constrained optimization problems using ODE systems. This approach was first introduced by [START_REF] Hopfield | neural" computation of decisions in optimization problems[END_REF] to solve the traveling salesman problem. Since then, neurodynamic optimization has been applied to a wide range of optimization problems, including linear and quadratic programming problems [START_REF] Xia | A recurrent neural network for solving linear projection equations[END_REF], general convex programming problems [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF], biconvex optimization problems [START_REF] Che | A two-timescale duplex neurodynamic approach to biconvex optimization[END_REF], global optimization problems [START_REF] Che | A collaborative neurodynamic approach to global and combinatorial optimization[END_REF], and stochastic optimization problems [START_REF] Tassouli | A neural network approach to solve geometric programs with joint probabilistic constraints[END_REF]. These methods typically use the Lyapunov stability theorem to prove that the constructed ODE system has a global convergence property. This means that any state solution of the ODE system converges to an optimal solution of the target problem.

In particular, neurodynamic optimization for solving pseudoconvex optimization problems has received widespread attention in recent years. Researchers have applied it to many applications, including portfolio optimization, energy efficiency optimization, and production planning [START_REF] Liu | A one-layer recurrent neural network for constrained pseudoconvex optimization and its application for dynamic portfolio optimization[END_REF][START_REF] Yang | Energy efficiency optimization in mimo interference channels: A successive pseudoconvex approximation approach[END_REF].

Various neurodynamic methods have been proposed to solve pseudoconvex problems with different types of constraints, such as linear equation constraints [START_REF] Guo | A one-layer recurrent neural network for pseudoconvex optimization subject to linear equality constraints[END_REF], bound constraints [START_REF] Liu | A one-layer recurrent neural network for constrained pseudoconvex optimization and its application for dynamic portfolio optimization[END_REF], convex inequality constraints [START_REF] Bian | Neural network for nonsmooth pseudoconvex optimization with general convex constraints[END_REF], and quasiconvex constraints [START_REF] Liu | A one-layer recurrent neural network for nonsmooth pseudoconvex optimization with quasiconvex inequality and affine equality constraints[END_REF].

Deep Learning for Solving Differential Equations

Another line of research included in our work is the use of deep learning to solve differential equations.

The idea of using neural networks to approximate the solutions of differential equations was first introduced by [START_REF] Dissanayake | Neural-network-based approximations for solving partial differential equations[END_REF], where training was performed by minimizing a loss function based on the network's satisfaction of the boundary conditions and the differential equations. [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF] showed that the network architecture could be designed to satisfy the boundary conditions, and this method was extended to systems with irregular boundaries [START_REF] Mcfall | Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions[END_REF].

With the advancement of deep learning, this approach has received renewed attention with the goal of solving high-dimensional nonlinear partial differential equations (PDEs) [START_REF] Han | Solving high-dimensional partial differential equations using deep learning[END_REF][START_REF] Huang | Partial differential equations meet deep neural networks: A survey[END_REF]. One notable approach is the use of physics-informed neural networks (PINNs) (Raissi et al., 2019a), which incorporate physical laws and boundary conditions into the network architecture and training process.

PINNs have been applied to a variety of engineering problems, such as fluid mechanics [START_REF] Samaniego | An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications[END_REF][START_REF] Cai | Physics-informed neural networks (pinns) for fluid mechanics: A review[END_REF]. Variations of PINNs have been developed to address different problem scenarios [START_REF] Lu | Physicsinformed neural networks with hard constraints for inverse design[END_REF][START_REF] Zhang | Learning in modal space: Solving timedependent stochastic pdes using physics-informed neural networks[END_REF][START_REF] Liao | Solving flows of dynamical systems by deep neural networks and a novel deep learning algorithm[END_REF] or to improve computational performance [START_REF] Jagtap | Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations[END_REF][START_REF] Yu | Gradient-enhanced physics-informed neural networks for forward and inverse pde problems[END_REF][START_REF] Sharma | Accelerated training of physics informed neural networks (pinns) using meshless discretizations[END_REF]. Software packages have been developed to facilitate the application of these methods (Lu et al., 2021a;[START_REF] Chen | Neurodiffeq: A python package for solving differential equations with neural networks[END_REF].

Contributions

Our main contributions in this work are as follows:

• We propose a novel approach that combines the advantages of neurodynamic optimization in modeling CPNOPs as ODE systems and the power of deep learning in approximating the solutions to these systems. By transforming the CPNOP into a neural network training problem, our solver eliminates the need for traditional optimization solvers or numerical integration methods.

• We design a specialized training algorithm that takes advantage of the problem structure of CPNOPs to optimize the performance of our proposed model. The neural network is trained to simultaneously satisfy the ODE system and minimize the CPNOP objective function.

• In our experimental results, we demonstrate that our proposed approach superiority over numerical integration methods in terms of solution quality and computational efficiency for solving CPNOPs. In addition, we show the performance of the proposed method in solving a large-scale CPNOP which is difficult to be solved by classical numerical integration methods.

Outline

The remainder of this paper is organized as follows: Section 2 provides essential background knowledge, including an introduction to CPNOPs and neurodynamic optimization. Section 3 presents our proposed neural network model and its application to optimization problems. The training process for the proposed neural network is described in Section 4. In Section 5, experimental results are presented and discussed, comparing our method to numerical integration techniques. Finally, Section 6 offers a summary and outlines future research directions.

Preliminaries

Constrained Pseudoconvex Nonsmooth Optimization Problem

We consider the following optimization problem:

                   min x f (x) s.t. x ∈ I = {x ∈ R n |g i (x) ≤ 0, i = 1, . . . , m}, x ∈ S = {x ∈ R n |Ax = b}, (1)
where x = (x 1 , x 2 , . . . , x n ) T ∈ R n are the decision variables, f : R n → R is the objective function, g(x) = (g 1 (x), g 2 (x), . . . , g m (x)) T : R n → R m is the function of the inequality constraints, and Ax = b is the equality constraints with A ∈ R p×n and b ∈ R p . n, m, and p denote the number of decision variables, inequality constraints, and equality constraints, respectively.

In this paper, we consider the case where f (x) is pseudoconvex and nonsmooth, g(x) is convex and smooth, and A is of full row rank. We denote x * as the optimal solution of the problem.

Assumption 1.

• There exists a point x ∈ R n such that x ∈ S ∩ int(I), where int(•) denotes the interior of the set.

• The objective function f (x) is Lipschitz continuous and regular on S ∩ I.

Neurodynamic Optimization

Now let y : R → R n be a time-dependent function. In this context, y(t) represents the state of the system at time t, and dy dt = Φ(y(t)) is an ODE system representing the dynamics of the system. A neurodynamic optimization method aims to construct an ODE system dy dt = Φ(y(t)) for the control of y(t) such that y(t)

settles down to an optimal solution of the considered CPNOP represented in equation ( 1).

Definition 1 (State solution). Consider an ODE system dy dt = Φ(y(t)), where Φ : R n → R n . Given an initial point y 0 ∈ R n , a vector value function y : R → R n is called a state solution, if it satisfies the ODE system dy dt = Φ(y(t)) and the initial condition y(0) = y 0 .

Definition 2 (Global convergence). Let y(t) be a state solution of an ODE system dy dt = Φ(y(t)), and let X * be the set of optimal solutions of the CPNOP (1).

The state solution of the ODE system is said to converge globally to X * , if y(t) satisfies

lim t→∞ dist (y(t), X * ) = 0,
where dist (y(t), X * ) = inf x * ∈X * |y(t) -x * |, and ∥•∥ is the Euclidean norm. If the set X * contains only one optimal solution x * , then lim t→∞ y(t) = x * , and the ODE system is globally asymptotically stable at x * .

In this paper, a modified version of the one-layer neurodynamic approach [START_REF] Liu | A one-layer recurrent neural network for nonsmooth pseudoconvex optimization with quasiconvex inequality and affine equality constraints[END_REF]) is used to model the CPNOP. We adapt the original approach to better fit the problem and improve its performance.

The resulting ODE system is given by:

dy dt ∈ -θ(t)(I -U ) m i=1 (1 -µ (g i (y (t)))) ∂f (y (t)) + ∂B(y (t)) -A T ρ(Ay (t) -b), (2)
where I ∈ R n×n is the identity matrix, and U = A T AA T -1 A. The components of this ODE system are described as follows:

θ : R → R is a function that controls the convergence of the state solution y(t) to satisfy the equality constraints, given by

θ(t) =      0, if t ≤ T 0 , 1, otherwise , (3)
where T 0 = 1+∥Ay 0 -b∥ 1 /λ min AA T , y 0 ∈ R n is an initial point of the ODE system, λ min AA T represents the smallest eigenvalue of the matrix AA T .

B : R n → R is defined as follows:

B(y(t)) = m i=1 max {0, g i (y(t))} (4)
Since every g i (•) is convex, B(•) is also convex, which means that ∂B(•) always exists. To compute ∂B(•), we first define an activation function µ : R → R as follows:

µ(s) =      1, if s > 0, 0, if s ≤ 0.
(5)

Then, the closed-form for ∂B(y(t)) can be given as follows:

∂B (y(t)) =                  0, if y(t) ∈ S ∩ int(I), i∈I 0 (y(t)) µ (g i (y(t))) ∇g i (y(t)), if y(t) ∈ S ∩ bd(I), i∈I 0 (y(t)) µ (g i (y(t))) ∇g i (y(t)) + i∈I + (y(t)) ∇g i (y(t)), if y(t) ∈ S\I, (6)
where I 0 (y(t)) = {i ∈ {1, 2, . . . , k}|g i (y(t)) = 0}, I + (y(t)) = {i ∈ {1, 2, . . . , k}|g i (y(t)) > 0}, and bd(•) denotes the boundary of the set.

ρ : R p → R p is an activation function defined as ρ(s) = (ρ (s 1 ) , ρ (s 2 ) , . . . , ρ (s p ))

T , where

ρ (s i ) =            1 if s i > 0, 0 if s i = 0, -1 if s i < 0. (7)
Theorem 1 [START_REF] Liu | A one-layer recurrent neural network for nonsmooth pseudoconvex optimization with quasiconvex inequality and affine equality constraints[END_REF]). Let Assumption 1 hold. Given any initial point y 0 ∈ R n , the corresponding state solution y(t) of the ODE system (2) globally converges to an optimal solution of the CPNOP (1).

Surrogate Model

Initial Value Problem Setup

In this work, we model the CPNOP in equation ( 1) using the ODE system outlined in equation ( 2). As stated in Theorem 1, given any initial point y 0 , the state solution y(t) approaches an optimal solution x * as time t goes to infinity, i.e., lim t→∞ y(t) = x * .

To determine the state solution, we must first specify an IVP by choosing an initial point and a time range. For simplicity, we set the initial point to be the all-zero vector, i.e., y 0 = 0. The time range is defined as the interval [0, T ], where T is a positive number given in advance.

We introduce the following definition:

Definition 3 (Final state). Given a state solution y(t) and a time range [0, T ], the state at the end of the time range, i.e., y(T ), is called the final state.

The final state y(T ) holds significance as it approximates the optimal solution x * , i.e.,

y(T ) ≈ x * . ( 8
)
The goal of this study is to effectively and accurately determine the final state of a given time range.

Neural Network Based Surrogate Model

We propose a neural network based surrogate model to approximate the solution of the constructed IVP.

The structure of the model is shown in Figure 1 (left). The equation below describes the surrogate model:

ŷ(t; w) = 1 -e -t N(t; w), t ∈ [0, T ], (9)
where N(t; w) is a fully connected neural network with trainable parameters w, and (1 -e -t ) is an auxiliary function that ensures that the initial condition ŷ(t = 0; w) = 0 is satisfied.

The auxiliary function (1 -e -t ) is a modified version of the Lagaris method [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF], where the exponential form has been shown to improve convergence in previous studies [START_REF] Mattheakis | Hamiltonian neural networks for solving equations of motion[END_REF]. In our experiments, we also observed similar benefits from using this function, likely because it reduces the influence of the neural network further from the initial time.

Approximate State Solution. The surrogate model is designed to approximate the state solution of the IVP constructed in Section 3.1. This is represented by the following relationship:

ŷ(t; w) ≈ y(t), t ∈ [0, T ]. ( 10
)
Predicted Soltuion to the CPNOP. In particular, we denote the final state of the surrogate model as:

x(w) = ŷ(t = T ; w), (11)
which serves as a prediction of the optimal solution of CPNOP. By combining equation ( 8) and equation ( 10)

at t = T , we have x(w) ≈ y(T ) ≈ x * . ( 12
)
Figure 1 (right) shows our proposed surrogate model as an approximate state solution, thus providing a prediction for the CPNOP. The CPNOP being considered in this figure has four variables. The colored lines represent the surrogate model, a 4-dimensional vector-valued function that approximates the state solution of the IVP. The red stars represent the final state, which is a 4-dimensional real vector that serves as the predicted optimal solution to the CPNOP.

Model Training

Loss Function

We define the loss function for the proposed surrogate model as follows: where Φ(•) refers to the ODE system outlined in equation ( 2), which corresponds to a CPNOP to be solved.

L(t, w) = ∂ ŷ(t; w) ∂t -Φ(ŷ(t; w)) , (13
∂ ŷ(t;w) ∂t
is the derivative of the model output ŷ(t; w) with respect to the input time t. This derivative can be calculated analytically using automatic differentiation tools such as PyTorch or JAX [START_REF] Paszke | PyTorch: An imperative style, high-performance deep learning library[END_REF][START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF].

In Figure 2, we show how the CPNOP is integrated into the loss computation, rather than being a component of the surrogate model structure. The surrogate model itself is an empty framework with no goal of solving a particular CPNOP. By reformulating CPNOP as an ODE system and embedding it in the loss function, we train the neural network to let the surrogate model satisfy the ODE system. This approach is similar to that used in PINN [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF], where physical rules are encoded as PDEs and incorporated into the loss function for training.

The goal of training the surrogate model is to minimize the following objective function: However, it is computationally intractable to compute E(w) due to the integral part. Therefore, in practice, we train the model by minimizing the following batch loss:

E(w) = T 0 L(t, w)dt, (14)
L(T, w) = 1 |T| t∈T L(t, w), ( 15
)
where T is a set of collocation time points, and |T| is the size of the set. In this study, we specifically adopt a uniform distribution over the interval [0, T ] for sampling collocation time points. While there is existing Algorithm 1 A deep learning-based solver for CPNOP Input: A CPNOP as defined in (1); A time range [0, T ] Output: Predicted optimal solution to the CPNOP 1: function CPNOP solver:

2:

Derive the ODE system, Φ(•), according to the CPNOP. ▷ By equation (2).

3:

Instantiate a surrogate model ŷ(t; w). ▷ By equation ( 9).

4:

Set ϵ best ← ϵ(x(w)) ▷ Initialize ϵ best .

5:

while iteration ≤ maximum iteration do 6:

Sample a batch of times T ∼ U (0, T ) uniformly from the interval [0, T ]. ▷ Data generation.

7:

Compute the batch loss L(T, w). ▷ Forward propagation.

8:

Update w by ∇ w L(T, w). ▷ Backward propagation.

9:

Project x(w) onto the feasible set I by x(w) ← P eq (x(w)). ▷ To satisfy equality constraints.

10:

Calculate the epsilon value ϵ temp ← ϵ(x(w)). 2023)), our practical experience indicates that the straightforward approach of uniform distribution yields comparable results.

Evaluation Metric

We propose a metric to evaluate the performance of the predicted optimal solution x(w). This metric, known as the epsilon value, is defined as follows:

ϵ(x(w)) =      f (x(w)), if x(w) ∈ S ∩ I, +∞, otherwise , (16)
where f (•), S, and I are given in equation ( 1). If x(w) is within the feasible set, the epsilon value is set to the objective value corresponding to x(w), i.e., f (x(w)). Otherwise, the epsilon value is set to +∞ to indicate that the prediction is not within the feasible set.

Since it is difficult to strictly satisfy the equality constraints, we use the following technique to project the prediction into S:

P eq (x(w)) = x(w) -A T AA T -1 (Ax(w) -b). ( 17
)

Pipeline

Algorithm 1 outlines the steps for solving a CPNOP using our proposed surrogate model approach. The algorithm starts by constructing an IVP from the given CPNOP and a pre-specified time range of [0, T ]. A neural network is then initialized as a surrogate model to serve as an approximate state solution for this IVP.

The parameters of the network are then optimized by performing gradient descent on the designed batch loss function, as defined in equation ( 15), in order to improve the accuracy of the approximation.

A key aspect of the proposed algorithm is the use of the evaluation metric, as defined in equation ( 16), to assess the performance of the model at each training iteration. The variables ϵ best and ϵ temp represent the lowest and current epsilon values achieved by the model, respectively. After each update of the network parameters, the algorithm compares ϵ temp with ϵ best . Suppose ϵ temp is less than ϵ best , indicating that the model found a better prediction. In this case, the algorithm updates ϵ best = ϵ temp and saves the current best prediction x(w). This design is similar to traditional optimization algorithms that check for optimality conditions (e.g., KKT conditions) at each iteration.

Experiments

Our proposed approach was evaluated on the Google Colab Pro+ platform, using PyTorch 1.12.1 with CUDA 11.2 for the neural network and JAX 0.4.1 [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF] for modeling the ODE system.

Section 5.1 presents a standard CPNOP, where we demonstrate the complete process of solving the problem using the proposed method and provide a detailed comparison between our approach and numerical integration methods. In Section 5.2, we showcase a large-scale CPNOP, which cannot be addressed by classical numerical integration techniques, and illustrate how our proposed method effectively tackles this challenge. Finally, Section 5.3 discusses the advantages and limitations of the proposed approach.

Comparing CPNOP Solutions: Our Approach vs. Numerical Integration

Example 1: Consider the following specific instance of a CPNOP:

min x f (x) = x 1 + x 2 + e |x2-1| -40 (x 1 + x 2 + x 3 ) 2 + 3 s.t. g 1 (x) = -3x 1 + 2x 2 -5 ≤ 0 g 2 (x) = x 2 1 + x 2 -3 ≤ 0 h(x) = x 1 + 2x 2 + x 3 -2 = 0 (18)
According to Theorem 7 in [START_REF] Liu | A one-layer recurrent neural network for constrained pseudoconvex optimization and its application for dynamic portfolio optimization[END_REF], the objective function in this instance is pseudoconvex. It is also non-smooth due to the presence of the absolute value function |x 2 -1|.

IVP Construction. We modeled this specific instance using the ODE system outlined in equation ( 2).

We set the initial point to y 0 = 0 and the time range to the interval [0, 10] to construct an IVP for this ODE system. The true state solution of this IVP is denoted as y(t), which is unknown and has no analytical solution. Our experiments showed that this particular instance of CPNOP results in a stiff ODE system, which leads to poor computational efficiency and accuracy when using traditional numerical integration methods. Figure 3: Comparison of solution processes between our approach and RK45. The top panel shows the evolution of ŷ(t; w) = (ŷ 1 (t; w), ŷ2 (t; w), ŷ3 (t; w)) T obtained by our approach at 0, 100, 1000, and 10000 TIs (from left to right). The bottom panel shows the evolution of ȳ(t) = (ȳ 1 (t), ȳ2 (t), ȳ3 (t)) T obtained by RK45 at 0, 1000, 5000, and 10000 CPs (from left to right).

Experimental Setup of Our Approach. To solve the instance presented in equation ( 18), we instantiated a surrogate model as follows: ŷ(t; w) = (1 -e -t )N(t; w), where t ∈ [0, 10], N(t; w) is a fully-connected neural network with one hidden layer of 100 neurons and the Tanh activation function. We then trained this surrogate model using Algorithm 1, with a maximum of 10,000 iterations, the Adam optimizer with a learning rate of 0.001, and a batch size of 128.

Experimental Setup of RK45. We performed a comparison of our proposed approach with the Runge-Kutta (RK) method, which is a commonly used technique among numerical integration methods. In particular, we chose the RK45 solver, which can be accessed via the Scipy library [START_REF] Virtanen | & SciPy 1.0 Contributors[END_REF], as a representative implementation of the RK method. We opted to present only a comparison with RK45 in Figure 3 and Table 1, as experimental observations indicate that the approximate state solutions obtained from different numerical integration methods do not show significant differences. We set the number of collocation points to 10,000, which were uniformly distributed over the time range of [0, 10]. The approximate state solution obtained by RK45 is denoted as ȳ(t). Both ŷ(t; w) and ȳ(t) aim to be as close as possible to the true state solution y(t) over the interval [0, 10].

Solution Process of Our Approach. In the top panel of Figure 3, we show the evolution of ŷ(t; w) as it approximates y(t). Our proposed approach achieves this by training the neural network with the parameters w. As the network parameters w are updated at each training iteration (TI), a new ŷ(t; w) is generated. The upper left corner of Figure 3 shows the initial approximation, while the upper right corner of the figure shows the final approximation after 10,000 training iterations, which is a more accurate approximation for the IVP and closer to y(t).

Solution Process of RK45. In the bottom panel of Figure 3, we show the evolution of ȳ(t) as it approaches y(t). RK45 achieves this by incrementally advancing the collocation point (CP) from the initial time t = 0 to the final time t = 10. The bottom panel of the figure illustrates how the approximation of ȳ(t) the quality of the solutions reported in Table 1, it is clear that our approach is able to produce faster and better solutions for the CPNOP instance in equation ( 18) than the numerical integration methods.

Large-Scale CPNOP

Example 2: Consider the following large-scale CPNOP:

min x f (x) = 999 i=1 x 2 i - x 2 i x i+1 s.t. g i (x) = x 2 10(i-1)+1 + x 2 10(i-1)+2 + . . . + x 2 10(i-1)+10 -20 ≤ 0, i = 1, 2, . . . , 100. h(x) = Ax -b = 0 (19)
This large-scale CPNOP is a modification from [START_REF] Bian | Neural network for nonsmooth pseudoconvex optimization with general convex constraints[END_REF], which has 1000 variables, 100 inequality nonlinear constraints, and 1 equality constraint. The objective function f (x) is pseudoconvex, and the inequality constraints are nonlinear and convex. In the equality constraint, b = 16, and the vector A can be accessed from the link1 .

IVP Construction. We construct an IVP by specifying the time range [0, 10] and an initial point y 0 , which can be accessed from the link2 . We use a fully connected neural network with four hidden layers, each containing 150 neurons, to solve this problem. Other training settings are consistent with Section 5.1.

Difficulty in Solving with Numerical Integration. The large-scale nature of Example 2 leads to an IVP with 1000 state variables, making it difficult to solve using classical numerical integration methods due to low efficiency. We attempted to solve the problem using methods such as RK45, but they took over an hour without providing effective predictions.

Constraint Violation Metrics. We define two metrics to describe the degree to which the predictions satisfy the equality and inequality constraints:

eq vio(x pred ) = ∥Ax pred -b∥ , ( 20
) ieq vio(x pred ) = g(x pred ) + , (21)
where g(x pred ) + = min{0, g(x pred )}.

Effectiveness of the Proposed Method in Solving Large-Scale CPNOP. Figure 4 shows the performance of our neural network method in solving Example 2. We can draw the following conclusions from the experimental results:

• As shown in Figure 4-(A), the loss function decreased from 175 to 20, indicating that the proposed model effectively solved the constructed IVP, even for large-scale optimization problems. The model provided an acceptable predicted state solution.

• The proposed method progressively approaches the optimal solution. As shown in Figure 4-(B), the objective value decreased from 17,000 to 2,500. Although the obtained predicted value is not the optimal solution under the current experimental settings, it is already very close and provides a meaningful upper bound for the original problem.

• The proposed method effectively finds feasible solutions to the problem. As shown in Figure 4-(C), after only 500 iterations, the equality constraint violation metric eq vio decreased from an initial 18,000 to 12, and the inequality constraint violation metric ieq vio decreased from an initial 15,000 to 0. Furthermore, the projection transformation (17) can adjust the predicted values to satisfy the equality constraints, reducing the eq vio metric to zero.

Discussion

Advantages. In comparison to numerical integration methods for solving CPNOPs, our approach presents three distinct computational benefits:

• Our method predicts the final state at each TI, whereas numerical integration techniques only offer a prediction upon completion of the entire solution process. This results in faster convergence and more efficient computation, as demonstrated in Tables 1 and2.

• By transforming the CPNOP into a neural network training problem, our approach allows us to capitalize on deep learning infrastructure and methodologies, enhancing computational performance.

• As evidenced in Section 5.2, our method is capable of solving large-scale CPNOPs, while classical numerical integration techniques demand significantly more computational time and still fail to deliver a feasible predictions Limitations. While our approach has demonstrated some advantages, it also has some limitations that must be acknowledged. One limitation is that it requires more complicated tuning of the hyperparameters to achieve optimal performance, unlike numerical integration methods, which are generally more straightforward and require minimal tuning. In addition, when dealing with a complex instance of CPNOP, such as those with many variables and constraints, it may be necessary to use a more advanced neural network architecture and allocate more resources to training in order to achieve a high-quality solution.

Large-Scale CPNOP. We believe that the proposed neural network approach may open new possibilities for solving large-scale CPNOPs. Traditional solvers, when applied to large-scale problems, often face inefficiencies due to memory limitations or iterative solving processes. In contrast, the method proposed here potentially circumvents these issues by transforming the problem into training a neural network. Training large-scale neural networks is a relatively common and well-researched task in the field of deep learning.

However, one challenge with the neural network method is the lack of guaranteed convergence. This is due to the fact that neural network training generally involves non-convex optimisation, where continuous loss reduction cannot be guaranteed. Nevertheless, we can use the neural network prediction as a starting point to warm start a traditional solver, thereby obtaining the convergent solution.

Conclusion

In this paper, we presented a novel approach for solving CPNOPs that combines neurodynamic optimization with deep learning techniques, specifically using deep learning-based differential equation solvers such as PINN. Our results demonstrated that this approach leads to faster convergence and superior solutions compared to traditional numerical integration methods. Our work establishes a link between CPNOPs and deep learning, opening up new possibilities in the field of nonlinear constrained programming.

However, it is important to note that the approach is still in its early stages and has limitations, such as potential reliability and robustness issues. In future work, we plan to improve the performance of our approach by exploring advanced deep learning-based PDE solvers and alternative neurodynamic optimization approaches, and to apply the method to other types of nonlinear programming problems.

Figure 1 :

 1 Figure 1: Left: The structure of the surrogate model. Right: The surrogate model as an approximate state solution on the time range [0, T ].

Figure 2 :

 2 Figure 2: Incorporating the CPNOP into the loss function

  which is an integral of the loss function over the pre-given time range [0, T ]. The loss value L(t, w) represents the error of the model at time t, and the objective function E(w) represents the total error of the model over the time range [0, T ].

Figure 4 :

 4 Figure 4: Performance of the proposed neural network method for Example 2: (A) Loss value, (B) Objective value, and (C) Constraint violation metrics. eq vio and ieq vio are defined in (20) and (21), respectively.

  PINNs exploring more sophisticated collocation point sampling methods to enhance computational efficiency (e.g., Tang et al. (2023); Nabian et al. (2021); Wu et al. (2023); Fang et al. (

	11:	if ϵ temp < ϵ best then
	12:	ϵ best ← ϵ temp .	▷ Update ϵ best .
	13:	xbest ← x(w).	▷ Update xbest .
	14:	end if
	15:	end while
	16:	return xbest
	17: end function
	research in

Table 2 :

 2 Evaluation of computational efficiency between our approach and various numerical integration methods. ms, s and h denote milliseconds, seconds and hours, respectively.

	Our approach	Numerical integration methods
	TI	Runtime CP	RK45 Runtime	RK23 Runtime	DOP853 Runtime	Radau Runtime	BDF Runtime	LSODA Runtime
	10	202ms	10	811ms	433ms	1.77s	1.06s	Fail	458ms
	100	893ms	100	980ms	1.23s	2.3s	1.33s	Fail	1.47s
	1000	8.47s	1000	11s	13.2s	3min 25s	>2h	Fail	20min 20s
	5000	40s	5000	3min 31s 5min 5s	29min 43s >2h	Fail	Fail
	10000 80s	10000 7min 45s 12min 15s 49min 55s >2h	Fail	Fail

https://drive.google.com/file/d/1duUyhCtW0qViVfX0VwBQcoiv9QTWvBeO/view?usp=share_link

https://drive.google.com/file/d/1V098TrlLgPH-WHrlgkX7njVLy-J9FBc2/view?usp=share_link

improves as the number of CPs increases, with the plot of ȳ(t) at 0, 1000, 5000, and 10,000 CPs corresponding to the solved time ranges [0, 0], [0, 1], [0,5], and [0, 10], respectively. 1: Evaluation of solution quality between our approach and RK45. The solution represents the predicted optimal solution for the CPNOP instance given in equation ( 18). The epsilon value serves as the evaluation metric, as defined in equation ( 16).

Comparison of Solution Quality. In Table 1, we present a comparison of the predicted optimal solutions for the CPNOP instance obtained using our approach and RK45, along with their corresponding epsilon values, where a lower epsilon indicates a better solution. The solution of our approach is represented by the final state of the surrogate model, i.e., ŷ(t = 10, w), where each TI corresponds to a different value of ŷ(t = 10, w). The solution of RK45 is given by ȳ(t = 10 * (CP/10000)), which is the final state of a solved time range. For example, when CP = 5000, the corresponding solved time range is [0, 5], and the solution of RK45 at this CP is the final state of this time range, i.e., ȳ(t = 5).

The results in Table 1 suggest the following:

• Our approach converges faster than RK45. In other words, our approach is able to find a satisfactory solution at an early stage of the solution process. For example, at the 100th TI, our solution has an epsilon value of -11.963, while the solution of RK45 only reaches an epsilon value of -7.872 at the 100th CP.

• Our approach yields a more accurate final solution than RK45. At the end of the solution process, i.e., at TI=10000 for our approach or CP=10000 for RK45, our solution reaches an epsilon value of -11.992, which is superior to that of RK45.

Comparison of Computational Efficiency. In Table 2, we evaluate the computational efficiency of our approach in comparison to a selection of popular numerical integration methods, such as RK45, RK23, DOP853, Radau, BDF, and LSODA, which are available through the Scipy library. The results are presented in terms of the CPU time required to solve the CPNOP, with lower values indicating superior efficiency. 2 indicate that our approach is significantly more efficient than the numerical integration methods. Specifically, our approach took only 80 seconds to complete 10,000 TIs, while the fastest numerical integration method (RK45) took more than 7 minutes to complete 10,000 collocation points. Some numerical integration methods, such as LSODA, and BDF, even failed to complete the task. Combined with

The results in Table


domain: chim sdv sdv.can

text:
Léna Carmès

Guillaume Bort

François Lux
email: francois.lux@univ-lyon1.fr

Léa Seban

Paul Rocchi

Zeinaf Muradova

Agnes Hagège

Laurence Heinrich-Balard

Frédéric Delolme

Virginie Gueguen-Chaignon

Charles Truillet

Stephanie Crowley

Elisa Belllo

Tristan Doussineau

Michael Dougan

Olivier Tillement

Jonathan D Schoenfeld

Needa Brown
email: ne.brown@northeastern.edu

Ross Berbeco
email: ross_berbeco@dfci.harvard.edu

AGuIX nanoparticle-nanobody bioconjugates to target immune checkpoint receptors

This article presents bioconjugates combining nanoparticles (AGuIX) with nanobodies (VHH) targeting Programmed Death Ligand 1 (PD-L1, A12 VHH) and Cluster of Differentiation 47 (CD47, A4 VHH) for active tumor targeting. AGuIX nanoparticles offer theranostic capabilities and an efficient biodistribution/pharmacokinetic profile (BD/PK), while VHH's reduced size (15 kDa) allows efficient tumor penetration. Site-selective sortagging and click chemistry were compared for bioconjugation. While both methods yielded bioconjugates with similar functionality, click chemistry demonstrated higher yield and could be used for the conjugation of various VHH. The specific targeting of AGuIX@VHH has been demonstrated in both in vitro and ex vivo settings, paving the way for combined targeted immunotherapies, radiotherapy, and cancer imaging.

Introduction

In recent decades, inorganic nanoparticles (NPs) have demonstrated significant promise for oncologic indications. One noteworthy example is the AGuIX (Activation and Guidance of Irradiation X) NPs. AGuIX are ultra-small NPs composed of gadolinium (Gd) chelates on a polysiloxane core, displaying a mean hydrodynamic diameter of 4 ± 2 nm. The corresponding mean molecular weight (MW) has been estimated at around 20 kDa. These NPs are currently being investigated in several clinical trials targeting different types of cancer, including brain metastasis (Phase II: NCT03818386 and NCT04899908), cervical cancer (Phase I: NCT03308604), glioblastoma (Phase I/II: NCT04881032), and pancreatic and lung cancers (Phase I/II: NCT04789486). 1,2 AGuIX NPs possess theranostic properties, providing dual functionality for imaging and radiotherapy. This is achieved through the presence of paramagnetic gadolinium ions (Gd3+) embedded within the NPs. These gadolinium ions contribute to positive contrast effects in T1 magnetic resonance imaging (MRI) and enhance the effectiveness of radiotherapy due to their high atomic number (Z=64). 3,4 The small size of AGuIX NP enables preferential passive uptake in tumors while facilitating rapid elimination through the kidneys(in humans t1/2 = 1.29 ± 0.27 h), thus minimizing the risk of toxicity. Additionally, AGuIX NPs demonstrate tumor uptake comparable to that of other organic or inorganic NPs, achieving approximately 1.19 ± 0.87% ID/g through passive targeting after 24 hours in large variety of animal model. 5,6 Nevertheless, to further enhance the efficacy of AGuIX NP, improving specificity and tumor retention time via active tumor targeting is an important next step for AGuIX translation and future generations.

Preclinical studies employing bioconjugation of AGuIX with antibodies and peptides have showcased the effectiveness of active targeting, resulting in a minimum 1.43-fold increase in tumor accumulation. 5,[7][8][9] Peptides have gained significant popularity as ligands for targeting cancer cells due to their compact size and ease of production. However, they often exhibit low affinity for the target compared to antibodies. 10 Conversely, antibodies have demonstrated remarkable efficacy as targeting molecules, but their larger size hampers tissue penetration and prevents access to certain areas in tumoral tissues. 11 For these reasons, our work focuses on bioconjugates based on nanobodies, otherwise known as Variable Heavy domain of Heavy chain (VHH). These targeting biomolecules have a molecular weight of approximatively 15 kDa and are roughly ten times smaller than monoclonal antibodies, and therefore offer enhanced efficiency in targeting. 12 Their reduced size enables them to readily penetrate tumor tissues and bind to a greater number of a receptors with enhanced affinity and specificity. Furthermore, VHH are swiftly eliminated from the bloodstream, minimizing the risk of toxic accumulation. They possess a remarkable folding capacity and robust physicochemical properties, endowing them with superior stability and excellent solubility. [12][13][14][15][16] These distinctive features enable VHH to overcome some limitations associated with monoclonal antibody therapies. 17,18 It is worth highlighting that both AGuIX NPs (nanoparticles) and VHHs (single-domain antibodies) exhibit a similar compact size, approximately 15 to 20 kDa, endowing them with similar in vivo properties, particularly in terms of tissue penetration, circulation, and clearance. By bioconjugating them, it is possible to combine these benefits in a single product. Herein, we design bioconjugates that target Programmed Death Ligand 1 (PDL1) and the Cluster of Differentiation 47 (CD47) receptors by using the A12 and A4 VHH respectively. [19][START_REF] Ingram | Proc. Natl. Acad. Sci[END_REF][START_REF] Ingram | [END_REF] A12 exhibits high affinity in the low nanomolar range for human and murine PD-L1, making it suitable for targeting its expression. Similarly, A4 exhibits a binding affinity of approximately 10 pM for murine CD47. [START_REF] Ingram | [END_REF][START_REF] Ingram | Proc. Natl. Acad. Sci[END_REF] The PD-L1 ligand and the innate immune regulator CD47 are relevant immunotherapy targets. PD-L1 can be expressed on tumor cells and binds to PD-1 on the surface of T-cells, triggering immune cell exhaustion and inhibiting anti-tumor immunity. Interfering with this negative immune checkpoint, PD-L1 inhibitors have demonstrated efficacy across many different cancer types. Specifically, atezolizumab, an anti-PD-L1 therapy, was approved for treatment of lung cancer in the U.S. in 2016 after showing a 25.4% overall survival rate after 3 years in the treatment of metastatic non-squamous or squamous non-small cell lung cancer (NSCLC). 22 Four years later, in 2020, duvarlumab, a second anti-PD-L1 antibody, was also approved for lung cancer treatment, demonstrating a 5-year overall survival rate of 42.9%, when administered following definitive chemoradiation therapy. 23 These two antibodies have been approved as a second treatment after one or two doses of chemotherapy. However, despite these promising results, the majority of patients with solid tumors do not respond to PD-L1 therapy, highlighting the need for combination approaches and better patient stratification.

24

CD47 is an important receptor expressed on tumor infiltrating macrophages that interacts with signal regulatory protein alpha (SIRPα). This regulator is overexpressed in various malignancies and prevents phagocytosis by providing a "don't eat me" signal to immune cells. Targeting and blocking CD47 on tumors restores this phagocytic response which can then upregulate secondary, adaptive responses.

20,25-27

Antibodies targeting the CD47 receptor are currently being evaluated in multiple clinical trials. [28][29][30] The limited efficacy of existing immunotherapy in treating most solid tumors has highlighted the need to better define the tumor immune microenvironment and develop combination treatment approaches. PD-L1 / CD47 AGuIX combinations could help achieve both aims. Three of the approved immune checkpoint inhibitors (ICI) that target the PD-1/PD-L1 axis require biomarker confirmation of PD-L1, highlighting the need for PD-L1 in vivo diagnostic tools. 31 Similarly, visualization of CD47 expression may help identify tumors amenable to macrophage based therapies. Both PD-L1 and CD47 inhibitors have been demonstrated to synergize with radiotherapy in preclinical models [32][33][34][35] ; enhancing the effectiveness of radiotherapy with AGuIX could further increase the potency of this combination. Therefore, our primary focus was on developing AGuIX NPs that could effectively target and delineate crucial immunomarkers, specifically PD-L1 and CD47. This research serves as a foundational demonstration of the chemical synthesis approach, aimed at creating a platform that harnesses both the diagnostic and therapeutic potential of AGuIX, in conjunction with the immunotherapeutic capabilities of immunocheckpoint targeting.

To achieve this objective, we grafted AGuIX NPs with two different VHHs for targeted applications using two highly specific techniques: 1) sortagging [36][37][38] and 2) click chemistry. [39][40][41] The sortagging is based on an enzyme specific reaction and click on an azide-alkyne reaction. 40,[START_REF] Baskin | Proc. Natl. Acad. Sci[END_REF][START_REF] Kim | [END_REF][44][START_REF]Explaining Nobel Prize-winning chemistry techniques[END_REF] A12 and A4-modified AGuIX were prepared to compare both conjugation techniques and to assess receptor targeting using various in vitro assays. For the selected coupling synthesis, the ex vivo approach was studied using A12 nanobody in a highly aggressive murine model of melanoma.

Materials and methods

AGuIX NPs

The Gd-chelated polysiloxane NPs (AGuIX) were provided by NH TherAguix (Meylan, France) as a lyophilized powder. Their synthesis has been extensively documented in the scientific literature. [START_REF] Le Duc | [END_REF]47 They contain roughly between 10 and 20 Gd atoms per particle, which can be quantified using ICP-MS (Inductively Coupled PlasmaMass Spectrometry). 3,48 The presence of the gadolinium atoms confers to AGuIX radiation dose amplification and MRI contrast properties. Throughout the paper, each AGuIX NPs concentration is stated in g L -1 of AGuIX NP or M of Gd element. Within the polysiloxane structure AGuIX contain primary amines in the (3-Aminopropyl)triethoxysilane (APTES) function. 49,50 Each NP is estimated to possess a similar number of amine functions as Gd chelates. 51,52 These amines will be utilized as the main functional groups for all the biofunctionalizations. The AGuIX NPs used in this study have all been pre-grafted on the free amino function present on polysiloxane matrix with the cyanine 5.5 fluorescent dye (detailed synthesis in ESI).

VHHs

The A12 and A4 VHHs were synthesized at the Massachusetts General Hospital (USA), as described a previously published protocol. AGuIX-Cy5.5-C(W/T)GGG (600 µM GGG(W/T)C final concentration) and VHH (30 or 40 µM final concentration) were introduced in buffer (50 mM Tris-HCl, pH 7.5, 150 mM NaCl). Ligation was triggered (t0) by the addition of StrA7m (2.5 µM final concentration). The mixture was stirred (orbital) at room temperature for 3 h. After reaction, the crude was negatively selected using Ni-NTA bead (400 µL slurry for 1 mg VHH, spin at 800 G for 1 min) to remove 6xHis-tag starting materials (StrA7m, VHH and released VHH residue after ligation), and filtered (MWCO 50 kDa vivaspin, dilution x100,000). AGuIX@VHH was purified by SEC (size exclusion chromatography, method B detailed in ESI) and stored at -20 °C. AGuIX-Cy5.5-C(W/T)GGG preparations and synthesis optimizations are detailed in ESI.

2.4. Synthesize AGuIX@VHH through click chemistry A12-azide (50 µM, final concentration in reaction) was added to AGuIX-DBCO at (100 AGuIX mg mL -1 final concentration in reaction) in 10 mM Phosphate-Buffered Saline (PBS) solution. The solution was stirred for 24 h at 4 °C. A temporal follow-up of the click reaction is monitored by SEC using method A (detailed in ESI). After reaction, AGuIX@VHH was purified by SEC (method B) to remove sub-products (NHS-terminated, azide-PEG4-NHS and AGuIX-DBCO excess) and stored at -20 °C. AGuIX-Cy5.5-DBCO and A12-azide preparations and synthesis optimizations are detailed in ESI.

Material characterizations 2.5.1. Size Exclusion Chromatography (SEC) UltraViolet/VisibleFluorescence

The Superdex 75 increase 10/300 GL (Cytiva) was used for SEC, and the process was conducted using an isocratic system of elution in acetate buffer 100 mM (pH 4.7) at a flow rate of 0.8 mL min -1 (25 μL injected). The G1311A pump (Agilent), Photodiode Array Detector G1315B (Agilent), and fluorescence detector G1321A (Agilent) were employed during the process to record absorbance at 295 nm and fluorescence at 280 nm excitation and 340 nm emission for the method A and fluorescence at 650 nm excitation and 665 nm emission for the method B. The data were acquired using on ChemStation B.04.03 SP1 and analyzed on GraphPad Prism 8.0.1. For preparative purification, the same column was employed within an AKTA PURE chromatography system (Cytiva) using an isocratic system of elution in phosphate-buffered saline (PBS) 100 mM at a flow rate of 0.8 mL min -1 (5 mL injected).

Induced Coupled Plasma Mass Spectroscopy (ICP-MS)

To achieve optimal results, the ICP-MS was operated under specific conditions, including a nebulizer gas flow of 0.84 L min -1 , plasma gas flow of 15 L min -1 , an auxiliary gas flow of 1.2 L min -1 , and a plasma radiofrequency power of 1600 W. The Syngistix 2.3 software was used to control the ICP-MS, and the tuning of all other parameters was carried out to optimize the Gd signal.

High-Performance Liquid Chromatography (HPLC)-ICP-MS

Flexar LC system (PerkinElmer) coupled with a Nexion 2000B (PerkinElmer) was used. The separation was executed using Superdex 75 increase 10/300 GL, and measurements were obtained through an isocratic mode of elution with acetate buffer 100 mM (pH 4.7) at a flow rate of 0.4 mL min -1 . The Gd signal was monitored utilizing isotopes 156 and 158 and Empower software version 7.3 was utilized to acquire the Gd signal.

Dynamic Light Scattering (DLS)

The hydrodynamic diameter distribution of the purified product was measured with a Zetasizer NanoS DLS instrument (laser He-Ne 633 nm) from Malvern Instruments.

2.5.5. Taylor Dispersion Analysis (TDA) TDA experiments were conducted using a TDA-ICP-MS hyphenation between a Sciex P/ACE MDQ instrument and a 7700 Agilent ICP-MS, described elsewhere. 53 Fused silica capillaries with an inner diameter of 75 µm and outer diameter of 375 µm, and a total length of 64 cm, were coated with hydroxypropylcellulose (HPC) using a solution of 0.05 g mL -1 in water. Detection was carried out by ICP-MS at m/z=158 with a data acquisition rate of 500 ms point-1. Between runs, the capillary was flushed at 5 psi for 5 min with the mobilization medium. Peak deconvolution was carried out using Origin 8.5 software. The detailed method is described in ESI.

Competition ELISA

High-affinity 96-well plates were coated with either 2 µg mL -1 PD-L1 or 1.25 µg mL-1 CD47 and allowed to adhere overnight at 4 °C. Plates were washed and blocked with a 10% FBS solution. Plates were incubated with sample VHH (A12 or A4) or NPs (AGuIX, AGuIX@A12, or AGuIX@A4) followed by incubation with 6.25 nM biotin-anti-PD-L1 or 50 nM biotin-anti-CD47. Plates were then incubated with avidin-HRP (Abcam ab7403, 1:40000 dilution) followed by TMB (ThermoFisher N301) to detect colorimetric changes. TMB conversation was stopped at 40 min and absorbance was read at 450 nm.

The equilibrium inhibitory dissociation constant (Ki) was curve fitted using Graphpad 8.0.1 Top and Bottom are the plateaus in the units of the y-axis. LogKi is the log of the molar equilibrium dissociation constant of the unlabeled ligand (AGuIX, AGuIX@A12, or AGuIX@A4). RadioligandNM is the concentration of the labeled ligand (biotin-anti-PD-L1 or biotin-anti-CD47) and HotKdNM is the equilibrium dissociation constant of the labeled ligand (biotin-anti-PD-L1 or biotin-anti-CD47).

Biacore analyses

Biomolecular interactions between immobilized receptor PD-L1 and CD47 and analyte AGuIX@A12 or AGuIX@A4 were assessed by Surface Plasmon Resonance (SPR) on Biacore 2000 instrument (Cytiva). 53 For the kinetic assays, the A12 was injected up to a concentration of 61.2 nM, the A4 up to 67 nM, the AGuIX NPs at 0-915.7 nM in Gd 3+ , the AGuIX@A12 VHH at 0-61.2 nM, (0-915.7 nM equivalent in Gd 3+ ) and the AGuIX@A4 VHH at 0-16.7 nM (0-240 nM equivalent in Gd 3+ ). The chip was prepared as described in SI. All characteristic interaction constants (equilibrium dissociation (KD) and kinetic rate of association and dissociation (kA and kD)) were determined by curve fitting using the Langmuir 1:1 binding model implemented in Biaevaluation software 4.1.1. 54 The detailed method is described in ESI.

Cell culture

Murine melanoma (wild-type B16F10 or hiPDL1-B16F10) cells were culture in 10% fetal bovine serum (FBS; Invitrogen, USA) and 1% pen/strep (10,000 U.mL-1 penicillin and 10,000 µg.mL -1 streptomycin; Invitrogen, US) supplemented Roswell Park Memorial Institute medium (RPMI 1640; Gibco, Invitrogen, USA) at 37°C, 5% CO2 and optimal humidity.

Fluorescence Microscopy

Cover slips were placed in 48-well plates and seeded with either hiPDL1-B16 or wild-type B16 cells (50,000 cells/well) and allowed to attached overnight at 37°C. Cells were incubated with Cy5.5 (683/703 nm) tagged NPs (0.2 mg.mL -1 ) for 1 h. Cells were fixed, blocked, and permeabilizing with a solution of 10% FBS and 0.3%

Ex vivo Biodistribution

The animals were anaesthetized with isoflurane (induction: 3%, maintenance: 1.5-2.0%) in a mixture of 100% O2 (flow rate = 1.0-1.5 L min -1 ). At day 8 post-inoculation, 1.92 ± 0.09 MBq (mean ± SD), corresponding to the dose of 7.12 ± 0.35 and 2.29 ± 0.07 µg of AGuIX (mean ± SD) for AGuIX-[ 89 Zr]-and for AGuIX-[ 89 Zr]@A12 respectively 55 , were intravenously injected (i.v.) into tumor-bearing mice (n = 8 for both AGuIX and AGuIX-A12, 22.0 ± 1.9 g per mouse). The mice were euthanized at 4 h and 24 h post-injection (n = 4 per time points) and tissue activity was determined for several organs of interest after harvesting (blood pool, intestines, kidneys, spleen, pancreas, liver, muscle, bone, brain, tumor). Activity in various organs of interest is represented in percentage of injected dose per gram (%ID/g) (Tab. S2).

Statistical analysis

The results are reported as mean ± SEM, as stated in the figure captions. Statistical analyses were conducted using GraphPad Prism (version 8.0.1). For comparisons involving three or more means, a one-way ANOVA followed by the Kruskal-Wallis nonparametric test (for non-Gaussian populations) was employed for in vitro internalization assays. All in vitro experiments were performed in triplicate. Mann-Whitney test for multiple comparisons was used to analyze ex vivo biodistribution data. Statistical significance was considered at p < 0.05.

Results and discussion

To evaluate the best approach for biofunctionalization of AGuIX NPs with VHH, we compared sortagging and click chemistry. The sortagging bioconjugation approach permits for the selection of specific sites for modification, enabling more accurate predictions regarding the biological impact of the chemical alteration. Modifying a protein terminus is expected to have a smaller effect on the protein's folding and functionality as well as the added advantage of defining the orientation of the immobilized protein, which is anticipated to better preserve its function compared to other immobilization techniques. This enzymatic reaction relies on the enzyme Sortase A (SrtA), a transpeptidase produced by Grampositive bacteria, which catalyzes the formation of a peptide bond between two peptides to attach specific proteins to the cell wall or pili assembly. [56][57][58] SrtA specifically recognizes the LPXTG sequence (with X representing any amino acid) and by nucleophilic attack makes it reactive to an N-terminal oligoglycine, forming a peptide bond between the threonine of LPXT and the oligoglycine. This natural ligation system has been repurposed for protein modification and has gained popularity as a research tool due to its high ligation selectivity, simplicity, robustness, and the availability of various SrtA variants (expressed in Escherichia coli or commercially available), as well as other required materials. 37,59 Sortagging reactions have been used to graft specific biomolecules, such as peptides or proteins, onto a protein or other biological targets of interest and has already proved efficacious in NP conjugation. 60-62 Click chemistry describes chemical reactions aimed at achieving high yield and high selectivity in the formation of carbon-heteroatom bond. Click chemistry reactions are based on the 1,3-dipolar cycloaddition of an azide and an alkyne to form a 1,2,3-triazole. This reaction has been widely utilized in various applications due to its simplicity in terms of preparation and purification steps, enabling rapid generation of new products with high reaction rates, such as bioconjugates. 60,[63][64][65] The formation of the triazole linkage is irreversible and quantitative, providing excellent reaction stability. Azide and dibenzocyclooctyne (DBCO) are selected as strainpromoted azide-alkyne cycloaddition reactions (SPAAC) tools due to their relatively small functional groups with a favorable rate constant (k2= 0.2-0.5 M -1 .s -1 )39 which have been shown to be effective. The crucial challenge of click chemistry is to preserve the integrity and functionality of the biomolecule after grafting. 39,66,67 Given the advantages and challenges associated with each method, we have synthesized AGuIX and VHH conjugates using both approaches to compare them.

Synthesis of AGuIX-VHH by sortagging

The sortagging reaction is based on enzymatic synthesis involving three primary steps: (1) introduction of maleimide (Mal) moieties on AGuIX NPs, (2) introduction of the peptide GGG(W/T)C using the Mal moieties, and (3) introduction of VHH through the transpeptidation enzymatic reaction (Fig. 1a). λem=340 nm) of the purified bioconjugate AGuIX@A12 (light red), AGuIX@GGG (grey) with A12 (dark red). (c) Measurement of hydrodynamic diameter by DLS of AGuIX@A12 (light red), AGuIX@GGG (grey) and AGuIX (black).

AGuIX NPs were modified to mimic the functionality of the Nterminal oligoglycine residue typically employed in sortagging. Initially, a Mal functional group was introduced onto the surface of AGuIX NPs using a bifunctional NHS/Mal linker. AGuIX-Mal was further modified with a short peptide linker that incorporates a cysteine amino acid at the C-terminal and a triglycine residue (GGG(W/T)C) at the Nterminal. To enable the use of sortagging, the VHH was engineered by incorporating a C-terminal StrA motif LPETGG and a 6xHis tag. The addition of the 6xHis tag motif facilitated protein purification and removal of the C-terminal residue released during sortagging. 20 The transpeptidation reaction between A12 and AGuIX-C(W/T)GGG was verified using SEC monitoring. After purification, the final product exhibited a retention time (tR) of 13.2 minutes, corresponding to a product size ranging from 29 to 44 kDa (Fig. 1b, Fig. S1.c). The equilibrium for the reaction was reached within 2-3 h at room temperature. Following isolation using Nibeads to capture and remove the 6xHis-containing reagents and by-products (StrA7m, residual peptides from VHH C-terminal, and unreacted VHH), AGuIX@A12 was purified via filtration and dispersed in PBS (Fig. S1.c). A VHH conversion of 32-33% was achieved during the sortagging process, resulting in an isolated yield of AGuIX@A12 of 5.3% (relative to A12) with a Gd/A12 ratio of 20 (Fig. 1b, Fig. S1.e). The hydrodynamic diameter (DH) of AGuIX@A12, as assessed by DLS, was 6.1 nm, slightly larger than that of AGuIX-C(W/T)GGG (4.1 nm) (Fig. 1c). This outcome confirms the modification of the particle surface and hints at the possibility of VHH grafting. After optimization, sortagging proved successful in grafting VHH onto the surface of AGuIX NPs. The key parameter to increase the grafting yield of VHH was to increase the equivalence of the oligoglycine in the reaction. However, this strategy also resulted in a higher NP/A12 ratio, which ideally should be around 1 for imaging applications. The NP/A12 ratio of 1 could be achieved after optimization of several parameters (Fig. S1.a to e), while VHH conversion decreased to 15%. The binding affinity, evaluated using the inhibitory constant (Ki) through a competitive enzyme-linked immunosorbent assay (ELISA), confirmed that the addition of AGuIX to A12 did not impact the competitive binding affinities (Fig. S2). Furthermore, AGuIX's lack of interaction with PD-L1 proteins also confirmed that the high affinity is solely attributed to the presence of the nanobody on AGuIX's surface. Moreover, the similar logKi values of 12.8 ± 1.1 nM, and 5.5 ± 0.6 nM for A12 and AGuIX@A12 obtained, respectively, strongly reinforced this result (Fig. S2).

Synthesis of AGuIX@VHH by click chemistry

Click chemistry based on strain-promoted azide-alkyne cycloaddition involves three primary steps: (1) introduction of the azide group on the VHH, (2) introduction of the DBCO group on the AGuIX NP, and

(3) introduction of the VHH through the click chemistry reaction (Fig. 2.a). The first step of the synthesis involves functionalizing the VHH with the azide group based on NHS-ester reaction with lysine residues. [68][69][START_REF] Nanda | Methods in Enzymology[END_REF] Functionalization was confirmed by SEC (Fig. S3.a) and matrix-assisted laser desorption/ionization-time of flight (MALDI-TOF) mass spectrometry (Fig. S3.b). 68 The AGuIX-DBCO was prepared by the same reaction consisting in NHS reaction with the primary amines on the NP surface. After synthesis and purification, the modified NP presented a ratio of 10 Gd/DBCO, i.e. 1 to 2 DBCO functions per NP (Fig. S3.c and d). AGuIX-DBCO high excess ratio (450Gd/A12) was chosen for the click chemistry to minimize any reactant rate limitations and ensure complete reaction (Fig. 2b). At least 95% of the VHH in solution successfully reacted with AGuIX, resulting in the formation of AGuIX@A12 with a size ranging from approximately 29 to 44 kDa, as determined by Superdex 75 protein calibration (Fig. 2d; Fig. S3.e and f). Purification by preparative SEC was effective in removing the unreacted AGuIXDBCO and other by-products, isolating the AGuIX@A12 bioconjugate. Dynamic light scattering (DLS) and Taylor Dispersion Analysis (TDA) size measurements of the bioconjugate were consistent (Fig. 3c, Fig. S3.i). The analyses of hydrodynamic diameters using DLS and TDA exhibit coherence and complementarity. DLS offers a standard resolution, providing an average value, but it is limited in its ability to differentiate between species as such ultrasmall sizes. On the other hand, application of TDA to nanoparticles is a recent method that has proven to be effective in distinguishing various populations of ultrasmall size. [START_REF] Labied | [END_REF] An increase in hydrodynamic diameter was observed following biofunctionalization (in DLS, DH, AGuIX = 3.0 ± 1.1 nm vs DH, AGuIX@A12 = 5.4 ± 3.1 nm). TDA is a highly accurate and absolute method based on deconvolution methods, enabling the hydrodynamic diameters of AGuIX NPs to be determined via their diffusion coefficients. 53 This analysis revealed the presence of two populations within the AGuIX@A12 sample in which 75 ± 2% of AGuIX NPs were effectively functionalized with the VHH, while 25 ± 2% remained unfunctionalized. The minority population (25 ± 2%) are the remaining AGuIX-DBCO (Fig. S3.j). In summary, AGuIX@A12 was obtained with a 20Gd/A12 (approx. 1 NP/A12) and a reaction yield close to 20% (based on VHH). Just as with the sortagging product, the affinity of AGuIX@A12 was assessed using ELISA, resulting in a logKi value of 13.0 ± 0.4 nM, similar to A12. The addition of AGuIX to the A12 nanobody through click chemistry did not influence competitive binding affinities; therefore, click chemistry did not compromise functional targeting properties (Fig. S2).

Methods comparison

First, we observed that both sortagging and click chemistry methods yielded AGuIX@VHH bioconjugates with relatively similar characteristics. The hydrodynamic size measurements (DH, sortagging: 6.1 ± 3.7 nm and DH, click: 5.4 ± 3.1 nm) and SEC chromatograms (tRsortagging: 13.2 min and tRclick: 12.9 min) were consistent between the two methods. Furthermore, a crucial factor determining the reliability of the synthesis is the maintenance of the VHH's affinity with the PD-L1 after its grafting to AGuIX. Evaluation of this affinity through a competitive ELISA demonstrated that, in both synthesis approaches, the binding affinity remained robust and similar to the reference VHH (logKisortagging: 5.5 ± 0.6 nM and logKiclick: 13.0 ± 0.4 nM). Although both conjugation products maintained similar binding affinities for the PD-L1 ligand, the reactivity of VHH upon click reaction conditions was higher, with a grafting yield before purification ηclick > 95% (vs ηsortagging < 50%). This observation aligns with existing literature, where sortagging reactivity efficacy may be limited due to equilibrium parameters (Fig S1 .e). 58,66,72 Moreover, the AGuIX-DBCO created during the click chemistry process could be separated from the final AGuIX@VHH bioconjugate, unlike the AGuIX-C(W/T)GGG used for sortagging. This separation provided a larger amount of purified AGuIX NP to be used for the click chemistry, reaching reaction equilibrium. Moreover, economic considerations are crucial when strategizing the scale-up of a process. It is noteworthy that biologics as used in the sortagging reaction are more expensive and less stable compared to chemical compounds. In conclusion, both sortagging and click chemistry methods afforded the AGuIX@VHH bioconjugates, however, the higher yield and lower cost associated with click chemistry have encouraged its adoption for further investigations (Table S1).

Proof of concept on A4 VHH

To demonstrate the reproducibility and robustness of the click chemistry approach, we conjugated a second VHH to AGuIX, this one specific to the CD47 receptor (A4. Mw =14.8 kDa). The first functionalization of the VHH with azide-PEG4-NHS demonstrated comparable efficacy to that of A12 (mean size of MWA12-azide: 15.33 kDa (2 azide grafted) and MWA4-azide:15.5 kDa (2-3 azide grafted)) (Fig. S6). The SEC indicate that the AGuIX@A4 (tR: 13 min) product has a size comparable to AGuIX@A12 (tR: 12.9 min) (Fig. 3). AGuIX@A4 was synthesized with a >95% conversion yield, similar to AGuIX@A12. After purification, the final purified product showed a ratio of 17Gd/VHH ratio, with a final yield of around 26%.

In vitro characterizations of AGuIX@VHH

A competitive ELISA confirmed that the addition of AGuIX to either VHH, A12 (Fig. 4a) or A4 (Fig. 4b), did not impact competition binding affinities. We demonstrated similar logKi values of 12.8 ± 1.1 nM and 13.0 ± 0.4 nM for A12 and AGuIX@A12, respectively, and 2.9 ± 0.2 nM and 3.3 ± 0.4 nM for A4 and AGuIX-A4, respectively. Moreover, AGuIX on its own exhibited no interaction with PD-L1 or CD47 ligands, providing additional support for the notion that the robust affinity observed primarily results from the incorporation of functional VHH on the AGuIX surface. These results were further validated using a second method: Biacore analysis which is based on a surface plasmon resonance (SPR) technique that which quantifies the association and dissociation phenomena. The affinity analyses conducted on the VHHs, and their associated receptors confirmed a strong binding affinity between these VHHs and their respective receptors, as shown by the equilibrium dissociation constants (KD) (0.63 ± 0.01 nM for A12 on PD-L1 and 0.70 ± 0.13 nM for A4 on CD47) (Fig. 4c, Fig. S4-5). These results

were supported by the literature, particularly in the case of A4, which has been extensively documented. 73 Comparing the KD values of AGuIX@VHH with their unmodified counterparts (0.39 ± 0.07 nM for AGuIX@A12 and 0.74 ± 0.01 nM for AGuIX@A4), we observed similar binding affinities further confirming that conjugation of AGuIX did not impact VHH binding. This similarity in KD values between AGuIX@VHH and unmodified VHH provided confirmation of the relevance of click chemistry in preserving the function of the VHHs. The A4, specific for the CD47 receptor, was selected to verify the specificity of the AGuIX@A12 to the PD-L1 receptor. The characteristic constants (equilibrium dissociation (KD) and kinetic rate of association and dissociation (kA and kD)) indicated a remarkable 10fold stronger and faster binding affinity for A12 towards PD-L1 compared to A4, regardless of the bioconjugation status. These results confirm the specificity of the interaction between the nanobodies and their respective receptors (detailed in ESI). To confirm ligand-receptor interaction on cells, internalization and receptor blocking assays on melanoma B16F10 tumor cells, modified for high PD-L1 expression, were conducted to further validate the functionality of VHH bound to AGuIX (Fig. 4.de). Melanoma hiPDL1-B16F10 murine cells are highly metastatic, aggressive models mimicking patient phenotypes, expressing both PD-L1 and CD47. 51 The conjugation of A12 and A4 significantly increased internalization of AGuIX as seen in the representative images (Fig. 4d) and corrected total cell fluorescence (CTCF, Fig. 4e). Blocking with either anti-PD-L1 or anti-CD47 decreased internalization confirming specificity and functionality of the AGuIX@VHH. Studies were repeated in wildtype B16F10 cells (lower PD-L1 expression) 74 (Fig. S7) with similar increase in internalization with the conjugation of A12 and A4 to AGuIX. Although complete loss of internalization was not observed, this may be an indication of only partial blocking of PD-L1 and CD47 receptors allowing for some receptor mediated internalization.

3.6.

Ex vivo AGuIX@A12 targeting Given the widespread use of PD-1 / PD-L1 inhibitors in clinical practice, we further evaluated the AGuIX@A12 product in proof-of-concept ex vivo studies. Targeting of PD-L1 was investigated in a murine melanoma model to confirm that AGuIX@A12 binding increasestumor accumulation. An ex vivo biodistribution study using zirconium-89 ( 89 Zr) radiolabeling was performed for accurate quantification of AGuIX and AGuIX@A12 in different organs at 4 h and 24 h post-intravenous injection (Fig. 4f -4i). A short and long time point was selected based on previous preclinical biodistribution studies and current clinical trials involving AGuIX. Currently, in clinical trials, radiotherapy is delivered 4 hours after AGuIX administration, justifying the clinical relevance of a 4-hour biodistribution time point and multiple preclinical trials showed that tumor accumulation was increased up to 4h, while the 24-hour time point was chosen to study the expected longer persistence of targeted nanoparticles in the tumor. 1,5,[75][76][77] As expected from previous work, both NPs accumulate within a few hours into the kidneys, leading to fast washout from the blood pool due to their ultrasmall size; lower accumulation of NPs in the liver can be detected. 76,78 Accumulation in the bone was attributed to the free [ 89 Zr] (non-coupled on the NPs) 79 . Both NPs can also accumulate in the tumor area due to the enhanced permeability and retention (EPR) effect, previously shown in several rodent models and in human studies. 6 However, a significantly higher accumulation of AGuIX-[ 89 Zr]@A12 was present in the tumor area compared to untargeted AGuIX-[ 89 Zr] (At 4 h, AGuIX: 2.10 ± 0.15 %ID.g -1 vs AGuIX@A12: 3.90 ± 0.67 %ID.g -1 ) (Fig. 4f -4i). This difference was observed 4 h post-injection and was more pronounced 24 h postinjection (at 24h, AGuIX: 1.34 ± 0.06 %ID.g -1 vs AGuIX@A12: 3.10 ± 0.15 %ID.g -1 ). This higher accumulation and retention of AGuIX@VHH suggests targeting of the PD-L1 receptors at the surface of cancer cells, that is expected to increase active accumulation and slow the AGuIX-VHH washout within the tumor area. This ex vivo validation of the higher targeting of AGuIX@A12 vs AGuIX is evidence of the in vivo targeting of PD-L1 by AGuIX@A12. AGuIX@A12 is thus a promising candidate for future work involving MRI-detection of PD-L1 or for more targeted radiosensitization of PD-L1 expressing tumors.

Conclusion

We report the preparation of AGuIX@VHH bioconjugates using two distinct methods, sortagging and click chemistry. The two methods resulted in NPs with similar physicochemical characteristics using the A12 VHH. However, click chemistry exhibited higher VHH conversion while using lower-cost reagents, making it a more favourable approach for further investigations. Interestingly, the two bioconjugation methods yielded AGuIX@VHH with strong binding affinities, suggesting that sitespecific VHH modification, such as sortagging, may not be required when coupling to relatively small NP such as AGuIX NP. PD-L1 targeting of AGuIX-A12 prepared by click chemistry was validated by ex vivo autoradiograpy demonstrating substantially greater uptake and persistence than untargeted NP. The A4 nanobody was used to the reliability and reproducibility of click synthesis for these constructs.

Overall, click chemistry emerges as a favourable and promising approach for preparing a broad array of potential AGuIX@VHH bioconjugates that can be applied to multiple nanobodies, with diverse applications in targeted therapeutic interventions and imaging. Further comprehensive evaluation and validation in larger preclinical studies and clinical trials are necessary to fully explore the imaging and therapeutic potential of AGuIX@VHH and their role in enhancing cancer treatment.

Fig. 1

 1 Fig. 1 Synthesis of AGuIX-VHH by sortagging reaction: method and characterizations (a) Scheme of the sortagging reaction, created with Biorender.com. (b) SEC chromatogram in the fluorescence intensity (λexc=280 nm;λem=340 nm) of the purified bioconjugate AGuIX@A12 (light red), AGuIX@GGG (grey) with A12 (dark red). (c) Measurement of hydrodynamic diameter by DLS of AGuIX@A12 (light red), AGuIX@GGG (grey) and AGuIX (black).

Fig. 2

 2 Fig. 2 Synthesis of AGuIX-VHH by click chemistry: method and characterizations (a) Scheme of the click chemistry reaction, created with Biorender.com. (b) SEC chromatogram in the fluorescence intensity (λexc=280 nm; λem=340 nm) of the purified bioconjugate AGuIX@A12 (red), AGuIX-DBCO (grey) and A12-azide (orange). (c) Comparison table of hydrodynamic diameters measured in TDA and DLS. Data presented as mean ± standard deviation.

Fig. 3

 3 Fig. 3 Reproducibility of the bioconjugation by click chemistry protocol on the bioconjugation of AGuIX with A4 VHH. The SEC chromatogram in the fluorescence intensity (λexc = 280 nm; λem = 340 nm) shows the purified bioconjugates AGuIX@A4 (light green), AGuIX-DBCO (grey) and A4-azide (blue).

Fig. 4

 4 Fig.4Binding affinity of AGuIX@VHH conjugates. ELISA plates were either coated with (a) PD-L1 or (b) CD47 proteins and affinity of free and AGuIX@VHH bound VHH was assessed. No difference in binding affinity for A12 or A4 free VHH compared to AGuIX validating that functional targeting capabilities of VHH are not impacted by conjugation. Control of AGuIX further confirmed that binding is due to presence of VHH conjugates rather than non-specific NP and protein interactions (n = 3). (c) Table of binding affinity constants for AGuIX@VHH conjugates. All characteristic interaction constants (KD, kA and kD) between the A12 and A4 to their respectively PDL1 and CD47 receptors were determined by Langmuir 1:1 binding model. Triplicates were made for all measurements except A4@AGuIX only one measurement was made. (d) Internalization of AGuIX@VHH in high PD-L1 expressing melanoma cells. Metastatic hiPD-L1 B16F10 cells were incubated with either Cy5.5 conjugated AGuIX or AGuIX@VHH (magenta) for 1 h prior to staining with CellMask plasma membrane (green) and DAPI (blue) nuclear stain. For a subset of samples, PD-L1 and CD47 receptors were partially blocked with commercial anti-PD-L1 and anti-CD47, respectively, prior to AGuIX@VHH incubation. (e) AGuIX and AGuIX@VHH fluorescence was quantified using ImageJ to calculate corrected total cell fluorescence, CTCF (n = 11 -40 cells, **** p < 0.001). (f,g) Quantitative analysis of ex vivo biodistribution (1.92 ± 0.09 MBq) at 4 h and 24 h after i.v. of (f) AGuIX@A12 and (g) AGuIX NPs which confirmed the specific active targeting of PD-L1 in B16F10 murine melanoma model. (n = 8; mean ± SD). (h) The tumor-to-blood and (i) the tumor-to-muscle ratios of AGuIX and AGuIX@A12 at different time points

  Zeiss AxioObserver microscope 63x. Fluorescence images were analyzed using ImageJ (version 2.14.0/1.54f) and corrected total cell fluorescence was calculated. [CTCF = Integrated DensityCell -(AreaCell x Mean fluorescence Background)] Glutamine, and maintained at 37 °C with 5% CO2 until confluence. Animal experiments have been performed according to the European directive 2010/63/EU and its transposition in the French law (Décret n° 2013-118). Experiments were conducted at the imaging facility CEA-SHFJ (authorization D91-471-105/ethics committee n°44). C57BL/6J mice (Janvier-Labs, France) were housed by 6 mice in each cage (bedding material: aspen wood) at room temperature 22°C, humidity 40%, under a regular 12-h dark/light cycle. Food and water were available ad libitum. 6 weeks old female C57BL/6 mice were purchased from Janvier laboratory. Mice were subcutaneously injected with 1 x 106 B16F10 cells suspended in DPBS (1× 100 µL) into both flanks while anesthetized with 2% isoflurane.

	2.10. Animal tumor model
	B16F10 cells (ATCC, USA) were cultured in DMEM media (Life Technologies, France) supplemented
	with 10% FBS (Life Technologies, France) and 1% penicillin/streptomycin (100 mg/mL), and 1% l-

Triton X-100. Cells were stained with CellMask Orange Plasma Membrane stain (ThermoFisher, 554/567 nm) and mounted with DAPI (350/470 nm) and Flouromount-G. Images were taken using a

Acknowledgements Special thanks to the international IRP Harvard-Lyon Radioboost project, which aims to boost radiotherapy with ultra-small nanoparticles. Support for this project was provided by Viewray, Inc (agreement #2019A017789).

Supporting Information

Supporting Information is available online or from the author.

Author Contributions

LC, GB and PR designed and performed the chemical experiments and product characterization. LS, ZM and NB conducted in vitro ELISA experiments and microscopy imaging. LH performed the Biacore experiments. AH conducted the TDA-ICP-MS experiments. VG-C carried out SEC preparative purifications. FD conducted the MALDITOF experiments. CT performed the ex vivo experiments. TD helped with NP formulation and characterization. LC, GB and NB wrote the manuscript. FL, OT, MD, JS, GB, NB and RB conceptualized the idea, supervised the study, and revised the manuscript. All authors discussed the results and approved the final version of the manuscript.

Conflicts of interest

LC is an employee of NH TherAGuIX that develops AGuIX NPs. OT and FL possess shares in this company. JS reports research support paid to the institution: Merck, BMS, Regeneron, Debiopharm, EMD, Seronoand consulting / Scientific Advisory Board / Travel fees: Castle Biosciences, Genentech, Immunitas, Debiopharm, BMS, LEK, Catenion, ACI Clinical, Astellas, Stimit, Merck KGA, SIRPant, EMD Serono. Stock options: Immunitas.