Joint Surrogate Learning of Objectives, Constraints, and Sensitivities for Efficient Multi-objective Optimization of Neural Dynamical Systems

Frithjof Gressmann; Ivan Georgiev Raikov; Seung Hyun Kim; Mattia Gazzola; Lawrence Rauchwerger; Ivan Soltesz

Figure 1: Surrogate-assisted optimization for multi-scale neural modeling at supercomputing scale. Top: overview of the iterative pipeline connecting biological data acquisition (A) , parameter sampling and simulation (B) , surrogate modeling (C) , multi-scale neural simulation (D) , and optimization/evaluation (E) . A , Experimental recordings are obtained from in vivo and in vitro preparations of biological systems. B , Reproducing these observations requires efficient exploration of a high-dimensional parameter space subject to multiple objectives (behavior matching) and constraints (experimental and biophysical). C , A transformer-based surrogate model maps parameter configurations to a shared latent representation and jointly predicts objectives (regression head) and constraint satisfaction (classification head), yielding a differentiable approximation of the search landscape. D , Multi-compartment neuron models with detailed ion-channel dynamics are assembled into morphology-preserving networks and evaluated in parallel on distributed high-performance computing infrastructure via MPI workers. E , The surrogate enables iterative model refinement by identifying regions of optimality across multiple objectives.

Figure 2: Neural network surrogate performance on single-cell optimization. A , To benchmark surrogate performance across diverse dynamical systems, 9 morphologically distinct CA1 interneuron types are jointly optimized over 10-15 parameters against 4 electrophysiological objectives (e.g., somatic firing response to current injection, objective #1 shown) subject to 7-8 constraints. Solid and dashed traces show optimized and un-optimized voltage responses, respectively. B , Standard space-filling sampling strategies (symmetric Latin hypercube, SLH; Latin hypercube, LH; Monte Carlo, MC; Sobol sequences) fail to approach the theoretical maximum normalized hypervolume, whereas Gaussian process regression (GPR)-based surrogate optimization consistently reaches it. C , All tested neural network surrogates (ResNet and FTTransformer variants with objective-only “o-” or joint constraint-and-objective “c+o”-heads) achieve Pareto-front

\varepsilon

-indicator values well within the 2% acceptable-solution threshold and outperform GPR and multitask exact Gaussian process (MEGP) baselines. D , This improvement is explained by lower absolute prediction error (normalized root-mean-square error, NRMSE) of the neural network architectures relative to the GP baselines throughout the optimization. E , ResNet variants achieve higher accuracy at lower inference cost than the GP baselines; the FTTransformer incurs moderately higher inference times but delivers the largest accuracy gains. Since accuracy gains translate to faster convergence, reduced evaluations relative to GPs outweigh any increase in inference overhead across all configurations. Statistics: in B , markers show the mean and error bars

\pm 1

SEM over

n=27

optimization runs per strategy (9 populations

\times

3 independent trials with distinct random seeds), and small dots the individual runs. In C , boxplot of the

\varepsilon

-indicator for

n=27

runs per method (9 populations

\times

3 trials). In D , boxplot across 81 fits (9 populations

\times

9 epoch checkpoints of one completed optimization run). In E , boxplot and dots show the individual timings of all

n=81

fits per model.

Figure 3: Joint differentiable representation of objectives, constraints and sensitivities. A , Left: schematic illustration how the joint surrogate gradient - combining objective and constraint partial derivatives - corrects the search direction away from infeasible regions toward optimal solutions. Right: normalized hypervolume over optimization epochs for representative in vivo and in vitro cell types. All surrogate-assisted strategies (neural network and GP baselines) converge faster and with fewer evaluations than full simulation alone (dashed line). Insets show normalized Inverted Generational Distance (IGD) to full-simulation solution at epoch 25; surrogate-assisted methods achieve

\leq 0.48\%

(in vivo) and

\leq 0.04\%

(in vitro) deviation from the full-simulation Pareto front. B , Global ranking of surrogate strategies across all single-cell problems (lower is better). Left: joint constraint-and-objective learning (c+o) on an FTTransformer backbone achieves the highest solution quality (IGD). Right: objective-only architectures (o-) converge fastest (hypervolume area under the curve, HV-AUC), revealing a trade-off between convergence speed and final solution quality. The FTTransformer backbone consistently outranks the GP baselines on both metrics. C , Surrogate-gradient-based sensitivity analysis (sgrad) yields IGD and normalized hypervolume comparable to established methods (Fourier Amplitude Sensitivity Test, FAST; Derivative-based Global Sensitivity Measure, DGSM) while being a direct by-product of the surrogate model, requiring no additional simulation budget. Statistics: in A , thick lines show the mean normalized hypervolume over

n=3

independent trials (distinct random seeds) and thin lines the individual trials; inset bars show the mean IGD

\pm 1

SEM with dots for the

n=3

trials. In B , bars show the mean rank

\pm 1

SEM over

n=10

optimization problems (9 CA1 populations and the motoneuron) and dots the per-problem ranks. In C , markers show the mean

\pm 1

SEM and dots the individual trials;

n=3

trials for IGD and

n=5

trials for normalized hypervolume.

Figure 4: Solving highly constrained problems with constraint-surrogate gradients. A , Left: characterization of constraint feasibility for a motoneuron optimization problem via random (Monte Carlo) sampling. While four of seven constraints are satisfied at appreciable rates (pre-spike count, 99.9%; resistance range, 94.8%; initial voltage, 94.8%; tau range, 40.8%), three constraints - first inter-spike interval (ISI), ISI adaptation, and monotonic frequency-current (F-I) relationship - yield 0% feasibility, preventing random sampling from discovering any valid solution. Right: computing the surrogate gradient with respect to the constraint classification head (

\nabla_{\mathbf{x}}f_{\mathbf{c}}

) guides the optimizer into feasible regions. Constraint-gradient-augmented strategies achieve rapid convergence in normalized hypervolume, outperforming both standard surrogate methods and the full-simulation baseline over 25 epochs. Random sampling (dashed black) fails to make meaningful progress. B , Left: analysis of the surrogate gradient-descent trajectory (blue) alongside its evaluation on the full simulation (green), confirming that the surrogate approximates the true objective landscape. Trace samples (orange dashed) are selected along the descent path by maximizing predicted diversity across all objectives. Right: augmenting the initial sampling set with 10 surrogate-descent trace samples substantially reduces Pareto-front prediction error (NRMSE) compared with the initial sampling alone or augmentation with 10 additional random samples, demonstrating the informativeness of gradient-guided trajectories for bootstrapping surrogate training. Statistics: in A (right), lines show the mean normalized hypervolume over

n=3

independent trials (distinct random seeds) and shaded bands

\pm 1

SEM. B (left) shows a single gradient-descent trajectory; B (right) reports the NRMSE of one surrogate fit per condition evaluated on 511 Pareto-front points, so no error bars are shown.

Figure 5: Joint-surrogate optimization of large-scale network problems. A , The multi-scale CA1 hippocampal network optimization problem spans biophysical scales (brain volume

\rightarrow

connectivity

\rightarrow

cell morphology) to computational scales (detailed micro-circuitry simulations), requiring high-fidelity, full-scale simulation. B , Total compute cost, measured in CPU-days (top axis) and wall-clock hours on 300 cores per worker (bottom axis), broken down into problem evaluations (solid bars) and surrogate training plus inference (hatched bars). The joint c+o-FTTransformer reaches solutions comparable to the MEGP baseline with 2x less compute and to the GPR baseline with 5x less compute. C , Normalized hypervolume over 50 optimization epochs. The c+o-FTTransformer (ours) and MEGP show similar initial convergence rates, both dominating the GPR baseline; after approximately 20 epochs the joint transformer advances to higher hypervolume, saving 54% of epochs relative to MEGP and 80% relative to GPR to reach equivalent solution quality. C-I Instantaneous population firing rates before (dashed) and after (solid) optimization for oriens-lacunosum moleculare (OLM) interneurons (dendritic inhibition, left axis) and pyramidal (PYR) cells (principal output, right axis). Optimization brings both populations into their respective target ranges (shaded bands), yielding balanced network activity. C-II Firing-rate deviation from target (0 = ideal) across four functional cell classes (principal output, somatic inhibition, dendritic inhibition, fast dendritic inhibition). The joint c+o-FTTransformer consistently achieves smaller deviations than the MEGP and GPR baselines, producing network dynamics closer to the biological target. All panels derive from a single optimization run per surrogate method (

n=1

) because of the computational cost. To gauge the run-to-run variability of the trajectory in C , we run

n=3

independent replicates (distinct random seeds) of the c+o-FTTransformer configuration for the first 5 epochs under a revised objective formulation whose absolute hypervolume values are not directly comparable to C , yielding two dimensionless spread proxies: (i) The SEM of the fraction of attainable hypervolume gained since epoch 0 is

0.011

at epoch 5, corresponding to

\pm 0.022

normalized hypervolume in C . (ii) The SEM of the hypervolume gain relative to its mean is

27\%

at epoch 5, corresponding to

\pm 0.001

at epoch 5 and

\pm 0.021

at the final epoch of C if the relative variability is held constant. Because the three methods in C start from the same initial sample, both proxies exclude initial-sampling variance and quantify the stochasticity of the surrogate optimization itself.

Condition	$q+k>100$	$q+k>200$
Max blocks $L$	2	2
Max embedding dim $d$	128	96
Max heads $H$	4	3
Min group size	2	4
Pool every $n$ blocks	1	1

Table S1: Adaptive scaling rules for FT-Transformer.

$w_{c}\!:\!w_{o}$	Constraint BCE $\downarrow$	Constraint AUROC $\uparrow$	Objective NRMSE $\downarrow$
0.25	$1.135\pm 0.078$	$0.993\pm 0.009$	$1.044\pm 0.156$
0.5	$1.060\pm 0.054$	$1.002\pm 0.011$	$0.960\pm 0.171$
1 (unweighted)	$1.000\pm 0.000$	$1.000\pm 0.000$	$1.000\pm 0.000$
2	$0.904\pm 0.085$	$1.011\pm 0.012$	$0.963\pm 0.161$
4	$0.822\pm 0.108$	$1.009\pm 0.023$	$0.873\pm 0.198$

Table S2: Effect of head weighting on the joint loss. Constraint and objective prediction under a 16-fold sweep of the constraint-to-objective loss ratio

w_{c}\!:\!w_{o}

, normalized to a constant sum. All values are ratios to the equally weighted model fit on the same data with the same seed with

1.0

denoting no change (mean

\pm

standard deviation over

n=9

fits from 3 CA1 populations

\times

3 seeds). The constraint cross-entropy (BCE) is reported at the converged validation value while AUROC and NRMSE are the best values scored on the set of candidates the optimizer proposed in the following epoch. The loss term the weighting targets moves monotonically across the sweep while neither accuracy does. Fits use the ResNet backbone at the production batch size on the accumulated samples (

N=2300

per population) of a completed run.

Category	Parameter	FT-Transformer	ResNet
Architecture	Number of blocks $L$ / $B$	3	2
	Block dimension $d$ / $d_{\text{block}}$	128	192
	Hidden multiplier $r_{\text{ff}}$ / $m$	2.0	2.0
	Number of heads $H$	4	–
	Embedding dim per head	32	–
	Pooling	CLS	–

Table S3: Surrogate hyperparameters. Default hyperparameters for both surrogate model architectures.

Package	Version	Purpose
Python	3.12	Runtime
dmosopt	0.78 (commit f034729 )	Optimization framework
MiV-Simulator	0.2.1 (commit e440eaa )	Network simulation
NEURON	8.2.6	Biophysical neuron simulation
TensorFlow	2.16.2	Neural network surrogate training
PyTorch	2.5.1 (CPU)	MEGP baseline (via GPyTorch)

Table S4: Software environment. Software dependencies and the exact versions used.

Parameter	Description	Lower	Upper
pp	Soma area fraction	0.1	0.9
Ltotal	Total length ( $\mu$ m)	10	80
gc	Coupling conductance (mS/cm 2 )	0.1	50
soma_gmax_Na	Soma Na + conductance (S/cm 2 )	0.001	0.9
soma_gmax_K	Soma K + conductance (S/cm 2 )	0.001	1.0
soma_g_pas	Soma leak conductance (S/cm 2 )	0	0.01

Table S5: CA1 interneuron parameter space. Parameter space for the BS (bistratified) interneuron model.

Population	$R_{\text{in}}$ range	$\tau$ range	$V_{\text{thr}}$	f-I steps	# params
Parvalbumin basket cells (PVBC)	[53, 179]	[5, 21]	$-37$	6	12
CCK basket cells (CCKBC)	[240, 320]	[22, 29]	$-43$	4	12
Ivy cells (IVY)	[167.5, 367.2]	[50, 70]	$-39$	4	12
OLM cells (OLM)	[500.5, 600.1]	[30, 50]	$-44$	4	12
Axo-axonic cells (AAC)	[65, 179]	[10, 14]	$-32$	5	12
Bistratified cells (BS)	[30.5, 109.1]	[11, 15]	$-32$	3	12

Table S6: CA1 electrophysiological targets. Electrophysiological targets for CA1 interneuron populations.

R_{\text{in}}

: input resistance (M

\Omega

);

\tau

: membrane time constant (ms);

V_{\text{thr}}

: spike threshold (mV). f-I column shows the number of current injection steps.

Parameter	Description	Lower	Upper	Lower	Upper
		Narrow range		Wide range
gc	Coupling conductance	0.1	2	0	2
soma_gmax_Na	Somatic Na +	0.1	0.3	0	1
soma_gmax_K	Somatic K +	0.01	0.3	0	1
soma_gmax_KCa	Somatic KCa	0.0001	0.01	0	1
soma_gmax_CaN	Somatic CaN	0.00001	0.03	0	1

Table S7: Motoneuron parameter space. Narrow range (used in the standard optimization) and wide range (used in Fig. 4 ).

Population	Slice count
Pyramidal (PYR)	7,825
Parvalbumin basket cells (PVBC)	142
CCK basket cells (CCKBC)	104
Ivy cells (IVY)	217
Neurogliaform cells (NGFC)	93
OLM cells (OLM)	43

Table S8: Neuron counts per population in the CA1 network slice used for surrogate optimization.

Experiment	Figures	Trials
CA1 single-cell surrogate optimization	2 , 3 , S2 , S3	3
CA1 no-surrogate baseline	3 , S3	3
CA1 optimizer comparison (NSGA-II, AGE-MOEA, SMPSO)	– (nadir estimation)	3
CA1 sensitivity analysis (surrogate gradient, FAST, DGSM)	3 , S5	5 (HV), 3 (IGD)
CA1 initial sampling strategy comparison	2 , S1	3
CA1 gradient-target comparison (BS)	S4	3

Table S9: Independent trials per experiment. Number of independent trials (replicates) per experiment. Each trial uses a distinct random seed generated by the experiment management framework ( machinable ).

Figure	Panel / plot	Centre, error, $n$
2 B	Sampling strategy comparison (HV)	Mean $\pm$ SEM; dots: individual runs; $n=27$ per strategy (9 populations $\times$ 3 trials)
2 C	Additive $\varepsilon$ -indicator	Bars mean $\pm$ SEM; dots: individual runs (runs beyond the axis drawn at its edge); $n=27$ per method
2 D	NRMSE / MAPE box plots	Median, IQR box, 1.5 IQR whiskers, dots individual fits within whiskers (fliers not drawn); $n=72$ – $77$ per model after outlier filtering of 81 (9 populations $\times$ 9 checkpoints, one run)
2 E	Inference time	Box plot as in D with individual timings; $n=81$ timed predictions per model
3 A	Convergence curves (HV vs. epoch)	Mean, thin lines individual trials; $n=3$ trials; inset IGD mean $\pm$ SEM with dots
3 B	Mean rank plot	Mean $\pm$ SEM over $n=10$ problems; dots: per-problem ranks

Table S10: Error representation across figures. Measure of centre, error representation and sample size

n

for every panel that reports summary statistics. Individual observations are overlaid as dots or thin lines wherever a mean is shown.

	gpr	megp	o	c+o	Mean
NGFC	$0.0112\pm 0.0076$	$0.0133\pm 0.0085$	$0.0081\pm 0.0021$	$0.044\pm 0.028$	$0.019\pm 0.012$
BS	$0.484\pm 0.055$	$0.553\pm 0.055$	$0.549\pm 0.082$	$0.520\pm 0.046$	$0.527\pm 0.059$
SCA	$0.85\pm 0.25$	$0.98\pm 0.14$	$0.81\pm 0.10$	$0.87\pm 0.36$	$0.88\pm 0.21$
OLM	$1.27\pm 0.14$	$1.22\pm 0.20$	$1.18\pm 0.19$	$1.23\pm 0.14$	$1.23\pm 0.17$
PVBC	$1.288\pm 0.045$	$1.266\pm 0.068$	$1.258\pm 0.052$	$1.284\pm 0.056$	$1.274\pm 0.055$
AAC	$1.364\pm 0.026$	$1.367\pm 0.033$	$1.365\pm 0.037$	$1.366\pm 0.026$	$1.366\pm 0.031$

Table S11: CA1 single-cell surrogate comparison. Normalized hypervolume (mean

\pm

std across

n=3

trials) for each surrogate strategy and cell type after initial sampling. GPR: Gaussian Process Regression; MEGP: multitask exact Gaussian process; o: objective-only neural networks (mean of o-ResNet and o-FTTransformer); c+o: joint constraint-and-objective neural networks (mean of c+o-ResNet and c+o-FTTransformer). Higher values indicate better Pareto front quality; the theoretical maximum is

1.1^{4}\approx 1.4641

. Populations are sorted by mean hypervolume (ascending). NGFC exhibits the lowest performance due to stringent constraint satisfaction requirements, while IS approaches the theoretical optimum across all methods.

Firing behaviour
		c+o-FTTransformer		c+o-ResNet
Constraint	Violations	AUROC	MCC	AUROC	MCC
Monotonic f-I	10704 (58.5%)	0.875	0.622	0.847	0.581
First ISI	11619 (63.5%)	0.881	0.652	0.855	0.616
ISI adaptation	11746 (64.2%)	0.875	0.638	0.853	0.610
Pre-stimulus spike count	314 (1.7%)	0.367	0.000	0.583	0.000

Table S12: Per-constraint feasibility prediction accuracy. Feasibility detection performance of the joint constraint head across the 9 CA1 interneuron populations for both backbones fitted on the same samples (computed within cells, then averaged). The joint surrogate is re-fit offline at nine epoch checkpoints of a completed optimization run and scored on the candidates the optimizer proposed in the following epoch. “Violations” are infeasible samples (and percentage) for that constraint across the

18{,}299

training samples. “all (conjunction)” scores the overall feasibility prediction. AUROC is threshold-free (chance =

0.5

) whereas for MCC

0

is no better than constant prediction.

Constraint pairs	Joint surrogate	One model per constraint
Within firing behaviour	$+0.733\pm 0.194$	$+0.892\pm 0.070$
Within measurement validity	$+0.894\pm 0.059$	$+0.698\pm 0.166$
Across the two groups	$-0.655\pm 0.132$	$-0.368\pm 0.195$

Table S13: Gradient geometry of the constraint predictions. Mean cosine similarity

\pm

standard deviation between the input gradients of the per-constraint feasibility losses, evaluated at the candidates the optimizer proposed and averaged over 3 random initializations

\times

9 populations within each group of constraints. Joint surrogate is the c+o-ResNet used throughout while one model per constraint fits a separate standalone constraint surrogate. Constraints are grouped as in Table S12 .

Model	IVY	PVBC	CCKBC	AAC	IS	SCA	BS	OLM	NGFC
No surrogate	$1.4624\pm 0.0007$	$1.4641\pm 0.0000$	$1.4637\pm 0.0002$	$1.4629\pm 0.0014$	$1.4636\pm 0.0001$	$1.4622\pm 0.0008$	$\mathbf{1.4506\pm 0.0037}$	$1.4640\pm 0.0001$	$1.4499\pm 0.0129$
GPR	$1.4618\pm 0.0006$	$1.4633\pm 0.0014$	$1.4637\pm 0.0001$	$1.4631\pm 0.0006$	$1.4635\pm 0.0001$	$1.4635\pm 0.0003$	$1.4461\pm 0.0057$	$1.4640\pm 0.0001$	$1.4171\pm 0.0078$
MEGP	$1.4618\pm 0.0014$	$1.4616\pm 0.0026$	$1.4639\pm 0.0001$	$1.4594\pm 0.0037$	$1.4633\pm 0.0001$	$1.4634\pm 0.0003$	$1.4115\pm 0.0283$	$1.4638\pm 0.0000$	$1.4459\pm 0.0234$
o -ResNet	$1.4623\pm 0.0014$	$1.4618\pm 0.0021$	$\mathbf{1.4640\pm 0.0000}$	$1.4639\pm 0.0001$	$1.4636\pm 0.0000$	$1.4636\pm 0.0001$	$1.4178\pm 0.0319$	$1.4640\pm 0.0000$	$1.4604\pm 0.0018$
c+o -ResNet	$1.4629\pm 0.0010$	$1.4632\pm 0.0015$	$1.4640\pm 0.0001$	$1.4639\pm 0.0002$	$1.4636\pm 0.0001$	$1.4636\pm 0.0001$	$1.3934\pm 0.0492$	$1.4640\pm 0.0001$	$1.4553\pm 0.0069$
o -FTTransformer	$\mathbf{1.4640\pm 0.0001}$	$1.4641\pm 0.0000$	$1.4639\pm 0.0000$	$1.4640\pm 0.0001$	$1.4635\pm 0.0001$	$1.4633\pm 0.0001$	$1.4158\pm 0.0375$	$\mathbf{1.4641\pm 0.0000}$	$1.4625\pm 0.0014$

Table S14: Final normalized hypervolume by population. Mean

\pm

standard deviation across three trials at epoch 25 (highest in bold).

Model	IVY	PVBC	CCKBC	AAC	IS	SCA	BS	OLM	NGFC
No surrogate	–	–	–	–	–	–	–	–	–
GPR	$19.1274\pm 7.9034$	$0.0496\pm 0.0097$	$0.5835\pm 0.3656$	$4.3584\pm 0.7938$	$0.6985\pm 0.1340$	$5.8518\pm 3.1411$	$\mathbf{0.1448\pm 0.0258}$	$0.7515\pm 0.2890$	$1.5082\pm 0.4431$
MEGP	$18.8594\pm 7.0845$	$0.1178\pm 0.0091$	$\mathbf{0.2477\pm 0.2582}$	$3.3730\pm 0.6306$	$4.0507\pm 0.3621$	$9.7238\pm 2.9401$	$0.3608\pm 0.0762$	$2.8427\pm 1.2167$	$6.1695\pm 1.1426$
o -ResNet	$9.4390\pm 3.0480$	$0.0661\pm 0.0067$	$0.7147\pm 0.2173$	$0.5974\pm 0.2301$	$\mathbf{0.1695\pm 0.0220}$	$3.7030\pm 2.0149$	$0.2370\pm 0.0257$	$1.3929\pm 0.5676$	$0.1316\pm 0.0031$
c+o -ResNet	$4.4141\pm 1.2012$	$4.0812\pm 0.0079$	$0.4014\pm 0.2736$	$1.9805\pm 0.9929$	$0.2327\pm 0.0064$	$5.2827\pm 2.1386$	$1.1835\pm 0.7002$	$1.5669\pm 0.6529$	$0.8262\pm 0.3534$
o -FTTransformer	$\mathbf{0.3824\pm 0.1547}$	$\mathbf{0.0007\pm 0.0004}$	$1.0022\pm 0.2273$	$1.9209\pm 0.6244$	$0.3261\pm 0.0349$	$0.1459\pm 0.1171$	$0.1820\pm 0.0258$	$\mathbf{0.7436\pm 0.4216}$	$\mathbf{0.1061\pm 0.0025}$

Table S15: Normalized inverted generational distance by population. Mean

\pm

standard deviation of the no-surrogate reference-front comparisons, normalized by the Euclidean norm of the population-specific nadir and reported as a percentage (lowest in bold).

Model	IVY	PVBC	CCKBC	AAC	IS	SCA	BS	OLM	NGFC
No surrogate	$35.053\pm 0.022$	$35.116\pm 0.022$	$35.095\pm 0.021$	$35.032\pm 0.060$	$35.124\pm 0.001$	$32.570\pm 1.628$	$\mathbf{34.435\pm 0.285}$	$34.616\pm 0.335$	$26.998\pm 1.232$
GPR	$35.069\pm 0.009$	$35.095\pm 0.040$	$35.107\pm 0.016$	$35.022\pm 0.084$	$35.123\pm 0.001$	$34.333\pm 0.792$	$34.107\pm 0.629$	$34.962\pm 0.105$	$33.183\pm 0.147$
MEGP	$35.050\pm 0.042$	$35.054\pm 0.052$	$\mathbf{35.121\pm 0.000}$	$34.972\pm 0.046$	$35.118\pm 0.003$	$34.463\pm 0.049$	$33.766\pm 0.653$	$34.651\pm 0.133$	$33.465\pm 0.923$
o -ResNet	$35.063\pm 0.023$	$35.063\pm 0.038$	$35.119\pm 0.013$	$35.088\pm 0.024$	$35.122\pm 0.001$	$34.735\pm 0.369$	$33.615\pm 0.865$	$34.972\pm 0.100$	$\mathbf{34.098\pm 0.066}$
c+o -ResNet	$35.088\pm 0.029$	$35.058\pm 0.037$	$35.118\pm 0.009$	$35.086\pm 0.002$	$35.122\pm 0.001$	$34.446\pm 0.324$	$33.354\pm 1.050$	$34.975\pm 0.105$	$33.679\pm 0.496$
o -FTTransformer	$\mathbf{35.089\pm 0.025}$	$35.116\pm 0.030$	$35.111\pm 0.015$	$35.093\pm 0.025$	$35.122\pm 0.002$	$\mathbf{34.827\pm 0.151}$	$33.771\pm 1.061$	$\mathbf{34.982\pm 0.107}$	$29.002\pm 5.353$

Table S16: Hypervolume area under the curve (HV-AUC) by population. Mean

\pm

standard deviation through epoch 25 (

n=3

). Higher values indicate faster convergence (highest in bold).

Figure S7: Constraint satisfaction under random sampling for the motoneuron benchmark. Heatmaps showing the percentage of random samples that satisfy each individual constraint (rows) and all constraints jointly (bottom row) for four sampling strategies (columns) in wide-range (left) and narrow-range (right) search-space configurations. Cell annotations give the percentage; color scale runs from 0% (dark) to 100% (light). Wide-range sampling yields near-zero joint feasibility, whereas narrow-range sampling achieves 0.02–0.05% joint satisfaction. Sampling strategy has a comparatively minor effect relative to search-space extent. Percentages are means over

n=3

independent sampling runs.

Figure S9: Surrogate gradient (sgrad) vs. black-box optimizer comparison. Heatmap of normalized hypervolume AUC (higher is better) for each combination of surrogate architecture (rows) and optimizer (columns) on the motoneuron benchmark over 10 epochs. Cell annotations report mean

\pm

s.d. across

n=3

independent trials. Color scale ranges from 0.95 to 1.0. “sgrad” denotes direct surrogate gradient optimization via backpropagation; “nsga2” denotes the standard NSGA-II black-box optimizer used in the main experiments. The sgrad approach achieves competitive AUC across surrogate architectures, demonstrating that for this benchmark, direct gradient-based candidate generation on the differentiable surrogate is a viable alternative to population-based black-box search. However, success of direct surrogate descent may be undermined if the solutions are found in regions of the search space the model predicts confidently wrong.

Model	IVY	PVBC	CCKBC	AAC	IS	SCA	BS	OLM	NGFC
MEGP (feasible only)	$1.4618\pm 0.0014$	$1.4616\pm 0.0026$	$\mathbf{1.4639\pm 0.0001}$	$1.4594\pm 0.0037$	$1.4633\pm 0.0001$	$1.4634\pm 0.0003$	$1.4115\pm 0.0283$	$1.4638\pm 0.0000$	$1.4459\pm 0.0234$
MEGP (all samples)	$\mathbf{1.4623\pm 0.0010}$	$1.4616\pm 0.0019$	$1.4639\pm 0.0001$	$\mathbf{1.4617\pm 0.0013}$	$\mathbf{1.4636\pm 0.0000}$	$\mathbf{1.4635\pm 0.0003}$	$\mathbf{1.4176\pm 0.0269}$	$\mathbf{1.4639\pm 0.0001}$	$\mathbf{1.4514\pm 0.0073}$
Rank-1 LMC	$1.4615\pm 0.0014$	$\mathbf{1.4630\pm 0.0009}$	$1.4635\pm 0.0002$	$1.4586\pm 0.0046$	$1.4635\pm 0.0001$	$1.4634\pm 0.0002$	$1.4144\pm 0.0211$	$1.4636\pm 0.0001$	–

Table S17: Final normalized hypervolume for the MEGP ablations. Values are mean

\pm

standard deviation across three trials at epoch 25. The highest mean in each population is shown in bold.

Model	IVY	PVBC	CCKBC	AAC	IS	SCA	BS	OLM	NGFC
MEGP (feasible only)	$22.1571\pm 8.4526$	$0.1339\pm 0.0362$	$\mathbf{0.3457\pm 0.3568}$	$\mathbf{4.0268\pm 1.1657}$	$3.2049\pm 2.3933$	$\mathbf{4.6331\pm 6.0008}$	$0.2789\pm 0.0894$	$\mathbf{2.4267\pm 0.7917}$	$1.9597\pm 3.0832$
MEGP (all samples)	$\mathbf{11.8538\pm 8.6729}$	$\mathbf{0.1183\pm 0.0518}$	$0.9549\pm 0.7954$	$4.0327\pm 0.6143$	$\mathbf{0.4153\pm 0.1145}$	$7.8556\pm 5.8246$	$\mathbf{0.1849\pm 0.0199}$	$4.4460\pm 3.2251$	$\mathbf{1.1435\pm 1.2977}$
Rank-1 LMC	$12.6766\pm 3.4653$	$0.1432\pm 0.0105$	$0.5509\pm 0.4716$	$4.6108\pm 0.9797$	$2.6264\pm 3.4252$	$8.0956\pm 2.0168$	$0.2512\pm 0.0764$	$4.4332\pm 3.1316$	–

Table S18: Normalized inverted generational distance for the MEGP ablations. Values are mean

\pm

standard deviation across three paired trial/reference-front comparisons, normalized by the Euclidean norm of the population-specific nadir and reported as a percentage. The lowest mean in each population is shown in bold.

Model	IVY	PVBC	CCKBC	AAC	IS	SCA	BS	OLM	NGFC
MEGP (feasible only)	$35.050\pm 0.042$	$35.054\pm 0.052$	$\mathbf{35.121\pm 0.000}$	$\mathbf{34.972\pm 0.046}$	$35.118\pm 0.003$	$34.463\pm 0.049$	$\mathbf{33.766\pm 0.653}$	$34.651\pm 0.133$	$33.465\pm 0.923$
MEGP (all samples)	$\mathbf{35.072\pm 0.014}$	$35.052\pm 0.048$	$35.106\pm 0.032$	$34.969\pm 0.048$	$\mathbf{35.124\pm 0.001}$	$34.297\pm 0.573$	$33.671\pm 0.750$	$\mathbf{34.826\pm 0.156}$	$\mathbf{33.822\pm 0.197}$
Rank-1 LMC	$35.054\pm 0.018$	$\mathbf{35.079\pm 0.029}$	$35.088\pm 0.032$	$34.971\pm 0.074$	$35.120\pm 0.003$	$\mathbf{34.801\pm 0.189}$	$33.585\pm 0.793$	$34.515\pm 0.088$	–

Table S19: Hypervolume area under the curve for the MEGP ablations. Values are mean

\pm

standard deviation across three trials through epoch 25. Higher values indicate faster convergence; the highest mean in each population is shown in bold.

Model	Fit (s)	Prediction (s)	Fit slowdown	Prediction slowdown
MEGP (feasible only)	$35.21\pm 29.36$	$0.35\pm 0.32$	–	–
MEGP (all samples)	$330.69\pm 68.53$	$2.00\pm 0.33$	$11.6\times$	$8.4\times$

Table S20: MEGP compute cost by training-set policy. Fit and prediction times are mean

\pm

standard deviation in seconds across the benchmarked population–epoch pairs. Slowdowns are medians of the paired ratios relative to feasible-only MEGP over the

n=9

population datasets (one fit per population at the final checkpoint of the Figure 2 inference analysis).

			All test points		Feasible test points
Epoch	Training set	$N_{\text{train}}$	NRMSE $\downarrow$	Spearman $\uparrow$	NRMSE $\downarrow$	Spearman $\uparrow$
5	MEGP (feasible only)	372	$0.44$ [0.32, 0.56]	$0.21$ [0.13, 0.27]	$0.26$ [0.21, 0.50]	$0.46$ [0.34, 0.54]
	MEGP (all samples)	1533	$0.36$ [0.20, 0.67]	$0.51$ [0.30, 0.55]	$1.95$ [1.44, 3.68]	$0.31$ [0.20, 0.44]
13	MEGP (feasible only)	607	$0.61$ [0.40, 1.15]	$0.24$ [0.23, 0.37]	$0.49$ [0.34, 1.20]	$0.55$ [0.45, 0.58]
	MEGP (all samples)	2333	$3.97$ [1.47, 9.90]	$0.29$ [0.28, 0.52]	$10.60$ [4.30, 82.92]	$0.05$ [-0.02, 0.28]

Table S21: Offline prediction accuracy of MEGP by training-set policy. At each checkpoint epoch the surrogate is re-fitted on all samples evaluated up to that epoch of the completed MEGP runs, either on the feasible subset only (the default) or on all samples including infeasible ones, and scored on the candidates the optimizer proposed in the following epoch. NRMSE (lower is better) and Spearman rank correlation (higher is better) are reported on all finite test points and on the feasible test points only; values are the median [interquartile range] over

n=9

population–trial pairs (9 populations

\times

1 trial).

N_{\text{train}}

is the mean number of training samples used by each ablation.

	Narrow range (MN)		Wide range (MN-r)
Method	HV-AUC $\uparrow$	Final HV $\uparrow$	HV-AUC $\uparrow$	Final HV $\uparrow$
Constrained BO ( $q=1$ )	$2.983\pm 1.114$ †	$1.4126\pm 0.0726$ †	$0.129\pm 0.026$ †	$0.0322\pm 0.0066$ †
Constrained BO ( $q=16$ )	$5.538\pm 7.831$	$0.7320\pm 1.0352$	$0.773\pm 0.159$	$0.0322\pm 0.0066$
Constrained BO ( $q=100$ )	$5.077\pm 5.468$	$0.9417\pm 0.8172$	$0.773\pm 0.159$	$0.0322\pm 0.0066$
GPR	$14.911\pm 0.351$	$1.4624\pm 0.0024$	$1.420\pm 0.141$	$0.0638\pm 0.0044$
MEGP	$15.332\pm 0.101$	$1.4640\pm 0.0000$	$17.431\pm 8.416$	$1.0236\pm 0.6091$

Table S22: Constrained Bayesian optimization vs. black-box strategies. Hypervolume area under the curve (HV-AUC; higher is faster convergence) and final normalized hypervolume. Values are mean

\pm

standard deviation across trials (best in bold). The

q=1

case runs a reduced budget by construction, so it is scored over its own shorter window (marked † ) and excluded from the bold comparison. For reference, the last two rows ( ‡ ) also list the constraint-gradient-augmented strategies of Figure 4 A, only run on the wide range configuration (

n=3

trials). Since they use the surrogate gradient, which none of the other methods do, they are also excluded from the bold comparison.

Joint Surrogate Learning of Objectives, Constraints, and Sensitivities for Efficient Multi-objective Optimization of Neural Dynamical Systems

Abstract

Figures & tables

Explore similar work

Constraint-Aware Discrete Black-Box Optimization Using Tensor Decomposition

DCEmbed: Scalable Optimization over Neural Surrogates

Expected Hypervolume Maximization for Multiobjective Optimization under Uncertainties