Antonio Lorente Mur
email: antonio.lorente-mur@creatis.insa-lyon.fr

Pierre Leclerc

Françoise Peyrin

AND Nicolas Ducros
email: nicolas.ducros@insa-lyon.fr

Single-Pixel Image Reconstruction from Experimental Data Using Neural Networks

Single-pixel cameras that measure image coefficients have various promising applications, in particular for hyper-spectral imaging. Here, we investigate deep neural networks that when fed with experimental data, can output high-quality images in real time. Assuming that the measurements are corrupted by mixed Poisson-Gaussian noise, we propose to map the raw data from the measurement domain to the image domain based on a Tikhonov regularization. This step can be implemented as the first layer of a deep neural network, followed by any architecture of layers that acts in the image domain. We also describe a framework for training the network in the presence of noise. In particular, our approach includes an estimation of the image intensity and experimental parameters, together with a normalization scheme that allows varying noise levels to be handled during training and testing. Finally, we present results from simulations and experimental acquisitions with varying noise levels. Our approach yields images with improved peak signal-to-noise ratios, even for noise levels that were foreseen during the training of the networks, which makes the approach particularly suitable to deal with experimental data. Furthermore, while this approach focuses on single-pixel imaging, it can be adapted for other computational optics problems.

Introduction

Single-pixel imaging is an extreme configuration of computational optics, where a single point detector is used to recover an image [START_REF] Edgar | Principles and prospects for single-pixel imaging[END_REF]. Since the seminal work by Duarte and coworkers [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF], single-pixel imaging has been successfully applied to fluorescence microscopy [START_REF] Studer | Compressive fluorescence microscopy for biological and hyperspectral imaging[END_REF], hyperspectral imaging [START_REF] Rousset | Time-resolved multispectral imaging based on an adaptive single-pixel camera[END_REF][START_REF] Arce | Compressive coded aperture spectral imaging: An introduction[END_REF], diffuse optical tomography [START_REF] Pian | Compressive hyperspectral time-resolved wide-field fluorescence lifetime imaging[END_REF], image-guided surgery [START_REF] Aguénounon | Single snapshot imaging of optical properties using a single-pixel camera: a simulation study[END_REF], short-wave infrared imaging [START_REF] Zhang | Dual-band single-pixel telescope[END_REF], and imaging through scattering media [START_REF] Li | Compressive ghost imaging through scattering media with deep learning[END_REF]. Single-pixel measurements can be modeled as dot products between an image and some two-dimensional functions that are implemented through a spatial light modulator [START_REF] Edgar | Principles and prospects for single-pixel imaging[END_REF]. To limit acquisition times, it is highly desirable to reduce the number of light patterns, which leads to an undetermined inverse problem to recover the image.

In the field of single-pixel imaging, deep learning has been used to unmix the fluorescence intensity and lifetime from time-resolved measurements [START_REF] Yao | Net-flics: fast quantitative wide-field fluorescence lifetime imaging with compressed sensing -a deep learning approach[END_REF][START_REF] Smith | Unmix-me: spectral and lifetime fluorescence unmixing via deep learning[END_REF]. Higham and coworkers [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF] proposed a convolutional auto-encoder for single-pixel image reconstruction imaging that outperformed compressed sensing approaches. This network directly maps the measurement vector to the desired image, using a fully connected layer followed by convolutional layers. Several Deep Learning architectures have been proposed since then, to solve single-pixel reconstruction problems. For instance, Li et al. in [START_REF] Li | Compressive ghost imaging through scattering media with deep learning[END_REF] used a similar network to that introduced by [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF]. For measurements with very low signal-to-noise ratios, in [START_REF] Rizvi | Deringing and denoising in extremely under-sampled fourier single pixel imaging[END_REF] the authors proposed a U-net for highly compressed single-pixel Fourier imaging, while in [START_REF] Hoshi | Single-pixel imaging using a recurrent neural network combined with convolutional layers[END_REF], a recurrent neural network was used. Further, Li and coworkers [START_REF] Li | Single-pixel compressive optical image hiding based on conditional generative adversarial network[END_REF] investigated a conditional generative adversarial network to reconstruct highly compressed data from measurements with random binary patterns.

In the present work, we propose a new deep learning methodology for image reconstruction from single-pixel measurements corrupted by a mixed Poisson-Gaussian noise model. Traditionally, inverse problems with data corrupted by Poisson noise are tackled using variance stabilizing transforms, such as the Anscombe transform [START_REF] Anscombe | The transformation of Poisson, binomial and negative-binomial data[END_REF], followed by a Wiener filter [START_REF] Wiener | Extrapolation, interpolation, and smoothing of stationary time series[END_REF]. However, the resulting image can be blurred, and in particular for undersampled data. More recent alternatives have been used to solve this issue by exploiting statistical or handmade image priors [START_REF] Lefkimmiatis | Poisson image reconstruction with hessian schatten-norm regularization[END_REF][START_REF] Ding | Statistical image reconstruction using mixed poisson-gaussian noise model for x-ray ct[END_REF][START_REF] Li | Deconvolution of poissonian images with the pure-let approach[END_REF].

Solving the resulting problems is usually prohibitive for real-time applications, as a single reconstructed image can take several seconds. In this regard, deep networks are ideal candidates due to their rapid evaluation. As such, deep learning has been used for single-pixel imaging in the presence of Poisson noise [START_REF] Liu | Photon-limited single-pixel imaging[END_REF], although no special modifications were made to the neural network to allow for the Poisson noise.

Many studies have explored many different neural network architectures; however, they usually consider convolutional layers that act in the image domain, while the raw data are in the measurement domain. Moreover, the interpretation of the output of a neural network is still an open question, specifically in the presence of noise where robustness issues can occur, as indicated by [START_REF] Kim | 3dsrnet: Video super-resolution using 3d convolutional neural networks[END_REF]. In particular, in the domain of single-pixel imaging, a neural network that is not finely tuned to the system can be outperformed by linear reconstructors [START_REF] Jiao | Does deep learning always outperform simple linear regression in optical imaging?[END_REF].

Contribution

First, we explore the best way to map the raw data to the image domain, before applying a cascade of convolutional layers. This mapping is traditionally learned or computed through the Moore-Penrose pseudo-inverse. We describe a linear mapping that can be interpreted in terms of traditional image reconstruction. This linear mapping is implemented as the first layer of a (nonlinear) deep neural network. We introduced this idea in [START_REF] Ducros | A completion network for reconstruction from compressed acquisition[END_REF] for noiseless data, and here we extend it to noisy data. In particular, we propose to estimate the image intensity, which is unknown in practice, and to use this estimation to approximate the covariance matrix of the measurements.

Next, we describe a machinery that allows a network to be trained in the presence of Poisson noise. We provide a normalization scheme that allows raw data with different orders of magnitudes to be considered, which is mandatory for realistic scenarios where the image intensity is not known.

Finally, we validate our approach by reconstruction of an experimental dataset that we acquired with varying integration times and light fluxes. Upon acceptance, the datasets will be made available alongside implementations of our reconstruction methods in the Python toolbox (SPyRiT [START_REF] Lorente Mur | Single-Pixel pYthon Image Reconstruction Toolbox (spyrit) Version 1.0[END_REF]).

Problem, assumptions and limitations

Single-pixel imaging

Let ∈ [0, 1] be the image to acquire. The main idea of compressive optics is to measure = 1 using hardware, and to recover using software. The system matrix 1 ∈ R × , with < , collects the patterns that are sequentially uploaded on a spatial light modulator, to obtain the measurement vector. The patterns are traditionally chosen from within a basis ∈ R × ; i.e., 1 = where = [ , 0] describes the subsampling strategy. Classical choices for the basis include Fourier, discrete cosines, wavelets, and Hadamard basis, as discussed in [START_REF] Ochoa | Assessing patterns for compressive fluorescence lifetime imaging[END_REF].

In the present study, we consider the Hadamard basis. The subsampling strategies can be divided into adaptive and nonadaptive. Adaptive strategies progressively determine during the acquisition [START_REF] Rousset | Adaptive basis scan by wavelet prediction for single-pixel imaging[END_REF], while nonadaptive chose beforehand. Here, for simplicity, we consider the nonadaptive strategy introduced in [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF][START_REF] Baldassarre | Learning-based compressive subsampling[END_REF].

Acquisition model

As Hadamard patterns contain negative values, we adopt a differential strategy [START_REF] Lorente Mur | Handling negative patterns for fast single-pixel lifetime imaging[END_REF]. We denote by m + the measurements of the positive parts of the patterns, and by m -the measurements of the negative parts. We model noise as a mixture of Poisson and Gaussian distributions [START_REF] Foi | Practical poissonian-gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Rosenberger | Emva 1288 camera characterisation and the influences of radiometric camera characteristics on geometric measurements[END_REF]. The Poisson noise is signal dependant and originates from the discrete nature of the electronic charge, while the signal-independent Gaussian noise accounts for readout and circuit fluctuations. For both the positive and negative measurements, we have m +,-∼ P ( +,-1

) + N ( dark , 2 dark ) (1)
where P and N are the Poisson and Gaussian distributions, is a constant that represents the overall system gain (in counts/electron), is the intensity (in photons) of the image (which is proportional to the integration time), dark is the dark current (in counts), and dark is the dark noise (in counts). We further hypothesize that dark and dark are independent of the image intensity .

The normalized measurements m are finally defined as

m = ( m + -m -)/( ), (2)
2.3. Deep learning for image reconstruction.

Deep-learning-based methods reconstruct an image ˜ using nonlinear mapping ˜ = G ( ) where the weights of the network are optimized with respect to a loss functions; e.g., to minimize the quadratic error over an image database

G = argmin 1 -1 =0 G ( ) - 2 . (3)
where ( ) 0≤ ≤ -1 are the image samples of the image database, and ( ) 0≤ ≤ -1 are the corresponding normalized noisy measurements given by Equation [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF]. Traditionally, neural networks are made of cascading layers, and can be written as

G = G • . . . • G 1 (4)
where G ℓ , 1 ≤ ℓ ≤ is the ℓ-th (nonlinear) layer of the network, and • is the function composition. The first layer is generally a fully connected layer that maps the normalized measurements ∈ R to a raw solution in * ∈ R . In Section 3.2, we explore different strategies for the design of this layer.

Theory

Experimental set-up

Our measurements are obtained from the single-pixel camera experimental set-up depicted in Fig. 1, as first described by [START_REF] Lorente Mur | Deep neural networks for single-pixel compressive video reconstruction[END_REF]. The telecentric lens (Edmund Optics 62901) is positioned such that its image side projects the image of the scene onto the digital micro-mirror device (DMD; vialux V-7001), which is positioned at the object side of the lens. The object is transparent and is illuminated by a LED lamp (Thorlabs LIUCWHA/M00441662). The DMD can implement different light patterns (denoted as 1 in Section 2) by reflection of the incident light onto a relay lens, which projects the light into an optical fiber (Thorlabs FT1500UMT 0.39NA). This optical fiber is connected to a compact spectrometer (BWTek exemplar BRC115P-V-ST1). For every object, we sequentially upload onto the DMD all of the = 4096 Hadamard patterns of dimension = 64 × 64 pixels. We can down-sample the full measurement vector a posteriori to achieve any sampling ratio. To consider different noise levels, we acquire the same object with varying integration times and neutral optical densities. . This set-up is composed of a sample (S) illuminated by a lamp (L) in front of a filter wheel (FW), a telecentric lens (TL), a digital micro-mirror device (DMD), some relay lenses (RL), an optical fiber (OF), and a spectrometer (SM).

Mapping of the raw data to the image domain

We propose to use a Tikhonov-regularized interpretable solution [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF] to map the raw data into the image domain. In particular, we choose

˜ = G 1 ( ˆ ) = * , (5)
where the regularized data * is such that

* = [ * 1 , * 2 ]
, where * 1 ∈ R and * 2 ∈ R ( -) are regularized versions of the acquired and missing coefficients, respectively.

Let be the variance of our noisy measurements, and are respectively the mean and variance of our prior model in the Hadamard domain. The computation of and is described in Section 3.3, and that of in Section 3.4. We derived the analytical Tikhonov regularized solution given by *

1 ( ) = 1 + 1 [ 1 + ] -1 ( -1 ), (6a) * 2 ( ) = 2 + 21 -1 1 [ * 1 ( ) -1 ]. (6b)
where -) are the blocks of the covariance in the measurement domain, 1 ∈ R and 2 ∈ R ( -) are the blocks of (see Section 3.3). Interestingly, Equation (6a) can be interpreted as the denoising of the raw data in the measurement domain, while Equation (6b) is the estimation of the missing coefficients from the denoised acquired coefficients.

1 ∈ R × , 21 ∈ R ( -)× and 2 ∈ R ( -)×(
To circumvent the difficulty of inverting the signal-dependant matrix in Equation (6a), we choose to neglect the nondiagonal terms of 1 ;Denoting 2 1 = diag ( 1 ), we get Fig. 2. Proposed network. The first two layers map the measurements into the image domain according to the analytical solution given by Equation ( 6). These two layers can be interpreted as a layer that denoises Equation (6a) of the raw measurements, followed by a layer that estimates the missing coefficients from the denoised measurements of Equation (6b). The raw image ˜ is corrected by a cascade of Image-domain convolution layers (CL) and non-linear layers, such as the ReLU layers.

Prior mean and covariance

We compute the mean and covariance as introduced in Equation ( 6) from the the samples of an image database. We consider the classical estimators of the mean and covariance

=       1 2       = 1 =1 , (8)
=       1 21 21 2       = 1 -1 =1 ( -) ( -) . (9)
Note that both quantities are computed once and for all. We load them into the graphical processor unit when the network is initialized, with no need to re-compute them later (e.g., during training or evaluation).

Noise covariance estimation

To use the mapping proposed in Section 3.2, we need an estimation of the covariance matrix of our measurements. Assuming that the raw measurements are independent, from Equation ( 2) we can determine the values of the covariance as

= Diag 2 = Diag 1 + 1 + 1 - 1 + 2 2 dark 2 2 , (10)
where Diag ( ) refers to the diagonal matrix where the diagonal coefficients are the elements of , and refers to the identity matrix. As 2 depends on the unknown image as well as on the intensity , we exploit the raw data that also depends on and . Recalling that the variance of a Poisson variable is the same as its expected value, we approximate the expected value by the noisy sample; i.e., 2 ≈ The experimental parameters dark , and dark can be estimated as described in [START_REF] Rosenberger | Emva 1288 camera characterisation and the influences of radiometric camera characteristics on geometric measurements[END_REF]. We describe how to estimate in Section 3.5.

1 2 ( ˆ + + ˆ --2 dark ) + 2 2 dark 2 2 . (11)

Estimation of the image intensity

For image denoising, estimation of can be achieved through fitting methods (e.g., see [START_REF] Foi | Practical poissonian-gaussian noise modeling and fitting for single-image raw-data[END_REF]). These methods usually exploit homogeneous regions of the image, and are therefore not suited to raw measurements. Instead, we propose a simple empirical method, which consists of estimating from a rough estimate of the nonnormalized image . Considering the pseudo inverse, which is a linear operator that scales with , we consider

˜ = max ∈ {1,..., } ( * ) , with * =       ˆ 0       (12)
In practice, this simple estimator can provide the order of magnitude of with good approximation, but it obviously leads to some inaccuracies. This is illustrated in Fig. 3, where we show the distribution of the error of our estimation over the STL-10 dataset [START_REF] Coates | An analysis of single-layer networks in unsupervised feature learning[END_REF] (i.e., 113,000 images obtained by merging the unlabeled, training and testing sets) for different values of (512, 1024) and (5, 50, 2500 photons). Overall, the relative error is usually below 50%.

Experiments

Training neural networks using simulated data

We train our network using the STL10 database [START_REF] Coates | An analysis of single-layer networks in unsupervised feature learning[END_REF], with = 105, 000 images that correspond to the 'unlabeled' and 'train' subsets. We consider 8,000 images for the testing (i.e., the 'test' subset). The original 96×96 images are resized to 64×64 using bicubic transform, and are normalized between -1 and 1.

We implement the full networks using Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] (version 1.5.1; cuda V10.2.89). We train the network by solving Equation (3) using the ADAM optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], with an initial learning rate of 10 -3 , which is halved every 10 epochs, for a maximum of 100 epochs. In our experiments, the validation loss traditionally stops decreasing after 70 epochs. The training phase takes 3 h and 45 min on a NVIDIA GP107GLM [Quadro P1000 Mobile]. As seen in Section 3.5, it is not realistic to consider the intensity as a known constant. Therefore, we make vary during the training, with mean and standard deviation . Therefore, the network has to learn how to be robust to vary about . To fit with the observations made in Fig. 3, we choose = 0.5 to account for the worst-case scenarios. As shown in Section 5, varying alpha during the training phase, i.e., considering different noise levels, is crucial to get a robust network.

However, such an approach does not guarantee that a network trained around a given intensity can perform well for another intensity. A neural network trained for a certain probability distribution cannot generalize 'for free' to another probability distribution. We will further develop this point in Section 5.

Our experiments were done with the same Image domain layer structure as in [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF]. This means that our network has three convolutional layers, with each layer separated by a ReLU layer. The first has a kernel size of 9 and a depth of 64, the second has a kernel size of 1 and a depth of 32, and the final one has a kernel size of 5 and a depth of 1. This structure can, however, be changed.

Experimental data

As shown in the first column of Fig. 4, we acquire three different objects: the LED lamp directly (first row); a cat from the STL-10 test set printed on a transparent sheet (middle row); and the Siemens Star resolution target (bottom row) [37]. The ground-truth image is computed (see Section 4.3) from a fully sampled measurement vector that is acquired with the highest signal-to-noise ratio (i.e., high flux illumination, long acquisition time). Specifically, we consider no neutral density, and the integration time to 1 ms per pattern for the lamp, to 4 ms per pattern for the STL-10 cat, and to 8 ms per pattern for the Siemens target,

For each object, we also acquire a high-noise dataset by placing neutral optical density behind the lamp, to reduce the light flux. We consider optical densities of 1.3 for the LED lamp, 0.6 for the STL-10 cat, and 0.3 for the Siemens target. We set the integration time to 4 ms for both the LED lamp and the STL-10 cat, and to 8 ms for the Siemens target, we retain only = 512 measurements for both the LED lamp and the STL-10 cat, which accelerates the acquisition by a factor of 8. For the Siemens target, which has a richer spatial frequency content, we keep more measurements ( = 2048; acceleration factor of 2).

Our estimated instrumental parameters were : dark = 1070 counts, = 1.54 count/electron, and dark = 53 counts.

Evaluation metrics

Given the ground-truth image , we compute the peak signal-to-noise ratio (PSNR) of a reconstructed image * as

PSNR( * , ) = 10 log 10 2 2 * -2 (13)
For the experimental data, we have no direct access to the ground truth image. This image is computed as = gt , where gt is a fully sampled measurement with high signal-to-noise ratio. Before computing the PSNR according to Equation ( 13), we normalize the ground-truth image in the range [-1, 1].

Results

Simulated results

Training for different levels of noise

In Table 1, we evaluate the effects of training a network under varying noise levels. We also simulate the acquisition of the test images with the same source intensity ( = 50 photons) and for varying intensities (mean values = 50 photons and standard deviation = 25 photons). Training a network under different noise levels is detrimental when the image intensity is known; on average, it results in a PSNR drop of 23.19 -22.13 = 1.06 dB.

On the contrary, noise-varying training improves the reconstruction when the image intensity is unknown; on average, we obtain an enhancement of 21.07 -19.33 = 1.74 dB (second row of Table 1).

This observation is crucial, as the exact image intensity is not available beforehand in real-life experiments. Although we can estimate the image intensity from the raw data (e.g., by the method introduced in Section 3.5), this inevitably leads to inaccuracies. In the following, we only consider realistic scenarios where the image intensity is not known, which requires training with varying noise levels. In Table 2, we report the reconstruction PSNRs obtained using the proposed network, and for a network with the same architecture where the mapping is learned ('Free Layer') for different levels of noise dictated by different image intensities ( = ∞ corresponds to the noiseless version of the proposed network). The standard deviation of the image intensity is set to 50%, in agreement with the findings of Section 3.5.

As expected, the network trained with no noise (i.e., = ∞) performs very poorly in the presence of noise. Therefore, we focus our analysis on the cases where our network is trained with noise. We divide our analysis into three cases, which depend on the relative noise levels used during training and testing.

Training noise and testing noise have the same levels

In Table 2, the PSNRs of these experiments are shown in blue (i.e., diagonals). In all of our simulations, the proposed method outperforms the 'Free layer' network proposed in [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF].

Training noise is higher than testing noise In Table 2, the PSNRs of these experiments are shown in red (below lower diagonal). We can observe that networks trained on high levels of noise (e.g., = 2 or 5 photons) with our proposed method perform poorly when they are tested on data with low levels of noise (e.g., = 50 or 2500 photons). For instance, a neural network trained with = 2 photons yields an average PSNR of 17.40 dB when it is tested for = 2500 photons. This is a PSNR drop of 4.77 dB compared to the optimal training conditions ( train = test = 2500 photons). We can also observe that in these cases the 'Free Layer' [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF] outperforms the proposed method. We can see that in the most extreme cases where the training was with very high levels of noise, and the testing was with very low levels of noise, the 'Free Layer' network can increase the PSNR by 20.09-17.40 = 2.69 dB on average.

Training noise is lower than testing noise

In Table 2, the PSNRs of these experiments are shown in green (above diagonals). We can see that the proposed network behaves similarly to the optimal training conditions. The worst drop in PSNR is relatively low, as 0.71 dB. This shows that the proposed network has high reconstruction quality provided that it is trained with relatively low levels of noise (e.g., = 50 or 2500 photons). Contrary to the previous scenario, the proposed method outperforms the 'Free Layer' network in all of these experiments. Here, the presence of the denoising layer has a great advantage; the proposed method generalizes better to high noise experiments.

Testing

Training 2. Reconstruction of peak signal-to-noise ratios (PSNRs) for different training strategies. 'Proposed' refers to the neural network method with the proposed mapping in Section 3. 'Free Layer' refers to the neural network method where the mapping of the raw measurements is learnt jointly with the postprocessing layers, as in [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF]. Note that the network trained with no noise ( = ∞) corresponds to the case where we choose = 0. From top to bottom, image acquisition is simulated assuming increasing light intensity (i.e., decreasing noise levels). From left to right, images are reconstructed by networks that are trained using decreasing noise levels. Blue font, the testing and training noise levels are the same; green font, the testing noise is lower than the training noise; and red font, the testing noise is higher than the training noise. To facilitate the comparison between the two networks, we underline similar PSNRs (i.e., difference <0.1 dB) and use bold font to indicate the best performing networks (i.e., difference >0.1 dB).

(in ph.) ∞ 2 ±

Experimental Data

Fig. 4 illustrates the performance of our networks on experimental data acquired with the experimental set-up of a single-pixel camera presented in Section 3.1. We compare the performance of our proposed method with the result of the total variation regularized solution [START_REF] Li | An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing[END_REF], the Tikhonov-regularized solution of Equation ( 6), the proposed method trained without noise ('Noiseless Net'), and with noise ('Proposed'), as well as the Tikhonov-regularized solution of Equation ( 6) to which we applied a state-of-the-art denoising method BM3D [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF] ('Tikhonov+bm3d'), and a neural network reconstructor with the same architecture as the proposed method but where the mapping is learned jointly with the postprocessing layer ('Free Layer'), as in [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF].

First, we observe that the network trained with no noise performs very poorly on noisy data, and it is outperformed by the pseudo-inverse in the case of the LED lamp and the STL-10 cat. Visually, the reconstructed image retains many of the artifacts introduced by Poisson noise, which emphasizes the need for accounting for noise when training a neural network. Furthermore, we observe that learning the fully connected layer leads to many artifacts and a loss of detail (see the star sector target with the 'Free Layer'). This can be attributed to the learning of the mapping for specific noise values, instead of adapting to the noise, as in our proposed method. We also observe that the Tikhonov-regularized reconstruction is itself quite powerful. Even if it keeps some of the reconstruction artifacts, it compares very well to the results obtained using total variation regularization. In terms of PSNRs, Tikhonov-regularized reconstruction performs very similarly to the 'Free Layer' neural network, which is strong motivation to use it for mapping to the image domain.

Finally, we observe that the proposed network and 'Tikhonov+bm3d' offer the smoothest reconstructed images that are almost free of the artifacts introduced by down-sampling or Poisson noise. In terms of PSNRs, we observe that our method outperforms the others on experimental data, while being similar to 'Tikhonov+bm3d'. However, the proposed method is 1,000 times faster than the latter (0.2 seconds for BM3D versus 0.2 milliseconds for the proposed network). Moreover, 'Tikhonov+bm3d' is a two-step method where denoising is agnostic to reconstruction, while our method solves both problems at once, adapting denoising to the noise level. Although the PSNR gain is relatively low compared to the Tikhonov-regularized reconstruction, the proposed network leads to images that are well enhanced visually.

Discussion

The proposed mapping generalises much better to higher levels of noise than the learned mapping proposed in [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF]. This is an advantage when dealing with experimental data in real time, as it can be time consuming to dynamically load different neural networks for different lighting conditions. In experimental scenarios, we do not know , and therefore our proposed mapping with a training value of = 2500 ± 1250 is the most suited for processing experimental data. Indeed, for all of the testing conditions, that neural network offers a very similar quality of reconstruction to the optimally trained neural network for that noise value.

Our method allows visually convincing reconstructions and good performance with respect to PSNR against several levels of noise, when trained for low levels of noise (high values of ). It therefore shows good proprieties for dealing with experimental data. Unlike previous single-pixel deep-learning studies, such as [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF][START_REF] Rizvi | Deringing and denoising in extremely under-sampled fourier single pixel imaging[END_REF][START_REF] Hoshi | Single-pixel imaging using a recurrent neural network combined with convolutional layers[END_REF][START_REF] Li | An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing[END_REF], we have showcased the robustness of our method to different lighting conditions, and given an interpretable meaning to the first layer of our neural network.

The noise level of the image under acquisition is a key feature in real-life experiments. This parameter can be estimated first, and then the network that fits the actual noise level can be evaluated. However, as the noise level is unknown, this requires the loading of several networks, which might be a severe limitation for real-time applications, and in particular if the models are too large to be all stored on the graphics processor unit. Our proposed method is robust to different levels of noise as long as it was trained with low levels of noise (high values of ). Therefore, the same network can be used almost optimally in a wide variety of situations. The Tikhonov-regularized mapping that we propose helps to denoise the raw data, to provide an approximate reconstruction to the convolutional layers that can focus on learning spatial features.

One limitation of our deep network compared to more classic approaches such as [START_REF] Lefkimmiatis | Poisson image reconstruction with hessian schatten-norm regularization[END_REF][START_REF] Li | Deconvolution of poissonian images with the pure-let approach[END_REF][START_REF] Luisier | Image denoising in mixed poisson-gaussian noise[END_REF]] is that there is no theoretical guarantee that it will work for any image; in particular, if the image under acquisition significantly differs from those of our training set. This is a common concern for deep-learning approaches that are, however, seen to work well in practice. While this tends to be confirmed by our experimental results where our approach worked on the Siemens Star sector and on a LED light (both of which are very different from the images of the stl-10 database), there are no theoretical guarantees that this will always be the case.

Another limiting aspect of our study might be the choice of our architecture, which is shallower than popular architectures, such as the U-Net used in [START_REF] Jin | Deep convolutional neural network for inverse problems in imaging[END_REF], and more recent variants. However, we are keen to keep the number of network parameters as low as possible, to keep both the training and evaluation times as short as possible. Another limitation of this study concerns the analysis of the PSNRs of the images reconstructed from the experimental data, where the ground truth is not known. We limit this common issue by acquiring fully sampled low-noise images. Finally, we only test our algorithms on × pixel images, with = 64. By considering a database with high-resolution images (e.g., ImageNet), our network can be generalized to handle the case where > 64, either directly or by implementing a patch-based strategy to limit the memory requirements.

Compared to most studies dedicated to deep learning for inverse problems, we mainly focus on the (fully connected) layers of the network, which acts in the measurement domain. The post-processing layers that act in the image domain can be replaced by any variants (e.g., U-Net, resNet, others). Moreover, our approach is compatible with architectures inspired by conventional variational methods (e.g., [START_REF] Aggarwal | Modl: Model-based deep learning architecture for inverse problems[END_REF]). This will be the object of future work. Although this work focuses on single-pixel imaging, it can be used for any linear problem where the measured data is scarcely sampled in an orthogonal basis, just by setting to the corresponding basis. Finally, this network can easily be adapted to other noise models, as we only need an estimation of the mean and covariance of the noisy measurements.

Conclusion

We present a deep learning method to recover an image in single-pixel imaging by introducing mapping of the raw data to the image domain. This mapping is implemented as the first layer of a neural network, which provides an approximate reconstruction to a traditional cascade of convolutional layers. We also describe a framework for training a network using simulated data only. We describe how to exploit all of the experimental parameters to deal with experimental data, and also how to normalize the measurement, which is fundamental when a nonlinear reconstructor (e.g., neural network) is considered.

Although our network is trained using simulations only, it performs very well on experimental data, even when the noise level is unknown. The estimation of the noise covariance coupled with an appropriate training process appears to be efficient in a wide variety of scenarios. This is particularly interesting in real-time configurations where it is not possible to evaluate many different reconstructors. In future work, we will explore more complex interpretable network architectures.

Fig. 1 .

 1 Fig.1. Optical set-up of the single-pixel camera[START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF]. This set-up is composed of a sample (S) illuminated by a lamp (L) in front of a filter wheel (FW), a telecentric lens (TL), a digital micro-mirror device (DMD), some relay lenses (RL), an optical fiber (OF), and a spectrometer (SM).

( 5 ,Fig. 3 .

 53 Fig.3. Box plot of the distribution of the relative error on the estimation of using Equation (12) with six combinations of values of ( , ).

Table 1 .

 1 Training and testing under varying noise levels. Top row: Data simulated for a given source intensity ( = 50 photons), which is the same for all of the test images.

	Testing		Training
	(in ph.)	50	50 ± 25
	50	23.19 ± 1.98 22.13 ± 1.88
	50 ± 25	19.33 ± 1.67 21.07 ± 1.63
	Bottom row: Data simulated for test images with varying intensities (mean	= 50
	photons, standard deviation	= 25 photons). Reconstruction PSNRs are reported for
	a network trained using constant intensity (middle column) and varying intensity (right
	column).

( ) = 1 +

1 /( 2 1 + 2 ) ( -1 ),(7)where division and multiplication apply element-wise.

Acknowledgments. This material is based upon work done on the PILoT facility (PILoT, INSA LYON, Bât. Blaise Pascal, 7 Av. Jean Capelle 69621 Villeurbanne).

Disclosures. The authors declare that they have no conflicts of interest.

Funding. This study was supported by the French National Research Agency (ANR), under Grant ANR-17-CE19-0003 (ARMONI Project), and performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon.

Data availability. Data underlying the results presented in this paper are available in Ref. [START_REF] Ducros | Single-Pixel Hyperspectral Imaging dataset Version 1.0[END_REF].

The following columns show reconstructions using the total variation regularized solution [START_REF] Li | An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing[END_REF], the Tikhonov-regularized solution of Equation ( 7), the noiseless proposed network (Noiseless Net) Section 3.2 (network trained with no noise and where we assume = 0), a Deep neural network with the same architecture as our proposed method, where the mapping is learned as in [START_REF] Higham | Deep learning for real-time single-pixel video[END_REF] (Free Layer), the Tikhonov-regularized method combined with BM3D [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF] denoising and the proposed network Section 3.2 (trained with = 50 photons and = 0.5 ). All of the PSNRs are computed as described in Section 4.3, with the first column as the ground-truth.


Fan Li

Wei Huang

Xingsheng Wang

Xiaojiang Lv

Fuhao Mo
email: fuhaomo@hnu.edu.cn

Effects of Active Muscle Forces on Driver's Lower Limb Injuries due to Emergency Brake in Various Frontal Impacts

Keywords: Lower limb, Muscle contraction, Frontal impact, Injury criteria, Finite element

Accident data shows that driver's kinematics response in real accidents can be significantly different from that in dummy or cadaver tests because of driver's muscle contraction. In this study, a finite element (FE) human body model consisting of an upper body of a dummy model and a lower limb-pelvis biomechanical model with 3D active muscles was developed to in-depth investigate lower limb injuries. Driver's emergency reaction during frontal impact was simulated by modelling muscle active contraction based on a series of volunteer experimental tests. Besides, a realistic impact environment with the response of the restraint system and the invasion of the driver's compartment was established in this study. The results show that muscle contraction can cause extra loads on lower limbs during the impact, which can increase the injury risk of lower limbs. As for the femur injury, muscle contraction caused an additional 1kN axial load on the femur and the femur resultant bending moment of active models was also higher by about 10-40Nm. Besides, the tibial index of the model with muscle activation was about 0.1 higher. Additionally, the results indicate that the femur injury is strongly related to the combined action of both axial force and bending moment. And the variation of the injury tolerance along the tibia shaft should be considered when evaluating the tibia injury. Overall, the current lower limb injury criteria can be still lack of robustness.

Introduction

With a wide use of restraint systems in cars, fatal injuries of head and chest have been greatly limited [START_REF]Traffic Management Bureau of Ministry of public security. Annual report of the People's Republic of China on road traffic accident statistics[END_REF]. However, lower limb injury prevalence has not been reduced [START_REF] Rudd | Updated analysis of lower extremity injury risk in frontal crashes in the United States[END_REF]. Previous biomechanical cadaver tests showed that most Knee-Thigh-Hip (KTH) injuries of lower limbs occurred in distal femur and patella [START_REF] Melvin | Impact response and tolerance of the lower extremities[END_REF][START_REF] Powell | Cadaver femur responses to longitudinal impacts[END_REF][START_REF] Viano | Mechanisms of femoral fracture[END_REF]. However, Rupp et al. [START_REF] Rupp | Biomechanics of hip injuries in frontal motor-vehicle crashes[END_REF] analysed the statistics of Japan traffic accident data and concluded that the pelvis joint injury accounted for 42% of the lower extremity injury in frontal crashes followed by the femur fracture accounting for 32%.

Other studies also showed similar results [START_REF] Kuppa | An overview of knee-thigh-hip injuries in frontal crashes in the United States[END_REF][START_REF] Rupp | The tolerance of the human hip to dynamic knee loading[END_REF]. To explain such difference, some studies indicated that previous biomechanical tests didn't properly consider the combined effects of femur axial forces and bending moments [START_REF] Laituri | Derivation and theoretical assessment of a set of biomechanics-based, AIS2+ risk equations for the knee-thigh-hip complex[END_REF][START_REF] Ivarsson | The tolerance of the femoral shaft in combined axial compression and bending loading[END_REF][START_REF] Untaroiu | A numerical investigation of mid-femoral injury tolerance in axial compression and bending loading[END_REF][START_REF] Untaroiu | A finite element model of the lower limb for simulating automotive impacts[END_REF]. Some believed that muscle active contraction during drivers' emergency braking action plays an important role [START_REF] Mo | Investigation of occupant lower extremity injures under various overlap frontal crashes[END_REF][START_REF] Chang | Development of a finite element model to study the effects of muscle forces on knee-thigh-hip injuries in frontal crashes[END_REF][START_REF] Chang | Predicting the effects of muscle activation on knee, thigh, and hip injuries in frontal crashes using a finite-element model with muscle forces from subject testing and musculoskeletal modeling[END_REF][START_REF] Mo | A Lower Limb-Pelvis Finite Element Model with 3D Active Muscles[END_REF].

Ejima et al. [START_REF] Ejima | Prediction of pre-impact occupant kinematic behavior based on the muscle activity during frontal collision[END_REF] drew a conclusion from Institute for Traffic Accident Research and Data Analysis (ITARDA) that more than 60% of car drivers in Japan made evasive manoeuvres on braking with their feet when they realized an unavoidable accident. In these emergency cases, they might prepare their bodies by bracing muscles in order to face the upcoming impacts. All these actions can change their kinematic responses during the crash due to the active muscular forces as well as the muscular stiffness. Therefore, occupant behaviours in real accidents can be significantly different from the responses obtained in dummy or cadaver tests.

Since no muscle active contractraction can be involved in the cadaver tests, it's an effective method to investigate active muscle effects in frontal crashes by using biomechanical human body models with muscle contraction modelled. Silvestri [START_REF] Silvestri | Development of a finite element model of the knee-thigh-hip of a 50th percentile male including ligaments and muscles[END_REF] developed a lower limb skeletal model for the study of lower limb injury mechanism based on the work of Schauer [START_REF] Schauer | Finite element modeling of the human anatomic pelvis and leg[END_REF]. This model used 1D spring elements and *MAT_SPRING_MUSCLE material in Ls-dyna to simulate muscles. Chang [START_REF] Chang | Predicting the effects of muscle activation on knee, thigh, and hip injuries in frontal crashes using a finite-element model with muscle forces from subject testing and musculoskeletal modeling[END_REF] estimated the activation levels of muscle groups by what is called equivalence principle based on electromyogram (EMG) signals. And the activation levels were then loaded into a lower limb FE model containing a 35-bundle 1D active muscle unit. However, muscles in his model were modelled with only 1D elements, the material parameters were antiquated, and the equivalence method to get the muscle activation levels needed to be verified.

With the development of modelling method, muscle models with higher fidelity gradually emerged. Chang et al. [START_REF] Chang | Development of a finite element model to study the effects of muscle forces on knee-thigh-hip injuries in frontal crashes[END_REF][START_REF] Chang | Predicting the effects of muscle activation on knee, thigh, and hip injuries in frontal crashes using a finite-element model with muscle forces from subject testing and musculoskeletal modeling[END_REF] developed a FE model of a midsize male's pelvis and his lower extremities with 35 Hill-type muscles in each extremity. He found that muscle activation increases the bending moment in the femoral shaft, thereby increasing the risk of femoral shaft fractures by 20%-40%. Both Behr [START_REF] Behr | Posture and muscular behaviour in emergency braking: An experimental approach[END_REF] and Iwamoto [START_REF] Iwamoto | Development of a human body finite element model with multiple muscles and their controller for estimating occupant motions and impact responses in frontal crash situations[END_REF] conducted volunteer emergency braking experiments and measured the surface EMG signal to obtain muscle activation levels, which were then used to simulate the lower limb muscle activity in fontal collisions. Generally, all of these models mentioned above used simplified loading conditions. The effect of the restraint system and the invasion of the driver's compartment were not taken into account. The specific injury criteria such as the femur axial force, the femur bending moment, or the tibial index were not given. And different frontal crash conditions were not considered.

Therefore, the present study aims to investigate the role of active muscles in the injury of occupant lower extremities during realistic frontal crashes. A lower limb-pelvis FE model with detailed 3D muscles was used in this study. Muscles' active levels during the emergency braking were obtained by volunteer emergency braking tests. Its biofidelity of both passive and active properties has been validated against recent experimental tests [START_REF] Mo | A Lower Limb-Pelvis Finite Element Model with 3D Active Muscles[END_REF]. A realistic frontal crash model for simulating real accidents was then established. Both the response of a detailed restraint system and the invasion of the driver's compartment were taken into account in this paper. Different kinds of frontal impacts with different overlap rates were simulated to investigate the effects of muscle active contraction during the emergency braking on lower limb injury occurrences.

Methods

Lower limb-pelvis FE model with active 3D muscles

A lower limb-pelvis model with active 3D muscles was used in this study to investigate the effect of muscle active contraction on lower limb injuries during emergency braking. The 𝐹 𝑐𝑒 = 𝐴(𝑡)𝐹 𝑙 (𝑙)𝐹 𝑉 (𝑉)𝐹 𝑚𝑎𝑥 [START_REF]Traffic Management Bureau of Ministry of public security. Annual report of the People's Republic of China on road traffic accident statistics[END_REF] Where 𝐹 𝑐𝑒 is the muscle active contraction force, A(𝑡) is the muscle activation level, 𝐹 𝑙 (𝑙) is the force-length curve, 𝐹 𝑣 (𝑣) is the force-velocity curve and 𝐹 𝑚𝑎𝑥 is the maximum isometric force. 𝐹 𝑙 (𝑙), 𝐹 𝑣 (𝑣) and 𝐹 𝑚𝑎𝑥 were given by empirical formulas summarized by relevant scholars. The specific definition of each parameter can be found in the previous research [START_REF] Mo | A Lower Limb-Pelvis Finite Element Model with 3D Active Muscles[END_REF]. To obtain driver's muscle activation level A(𝑡) during the emergency braking, volunteer emergency braking experiments were conducted in Laboratory of Applied Biomechanics of IFSTTAR. The volunteers were asked to sit on the simulator seat according to the driver's sitting posture, and keep a constant 4N force on the acceleration pedal, simulating an average driving speed of 120km/h. In order to simulate the emergency braking situation, after receiving the braking order, the volunteers slammed on the brakes for at least 3 seconds.

Braking force, braking time and muscle EMG activities were recorded during the experiments.

Details regarding the experimental procedure and experimental data processing have been reported in the previous studies [START_REF] Mo | A Lower Limb-Pelvis Finite Element Model with 3D Active Muscles[END_REF][START_REF] Jammes | Emergency braking is affected by the use of cruise control[END_REF]. The average muscle active levels obtained after EMG signal data processing are shown in Table 1. 1. More details supporting the present model have been reported in the previous study [START_REF] Mo | A Lower Limb-Pelvis Finite Element Model with 3D Active Muscles[END_REF]. Then, a realistic frontal impact model including a detailed restraint system was eventually established based on the aforementioned driver's emergency reaction as shown in Fig. 2. Three typical frontal crash situations with different overlap rates were simulated including 100% full impact, 40% offset impact and 25% offset impact. As shown in Table 2, in order to investigate the effect of active muscles on lower limb injury, the simulations were classified into two groups with or without muscle activation. Each simulation consumed 24 CPUs and 8 hours. If the hybrid dummy and a full vehicle were both included to complete the whole collision process, the calculation process would be very long. Therefore, in this study, the acceleration curve of the B pillar and the invasion curve of the driver's compartment under different collision conditions were used as the input conditions of the simulations.

Establishment of realistic impact environment

100% frontal rigid wall collision (50km/h) and 40% offset collision (64km/h) based on China-New Car Assessment Program (C-NCAP) and 25% small offset collision (64km/h) based on The US Insurance Institute for Highway Safety (IIHS) were simulated with a full vehicle model. Under different frontal collisions, the acceleration curves of B pillar (Fig. 3)

and the invasion curves of the driver's compartment (Fig. 4 and Fig. 5) were recorded as input conditions of simulations in Table 2. Taking the 100% full impact simulation as an example, muscle active process and the input condition were set as shown in Fig. 6. When drivers are aware of the danger, most of them begin to slam on the brake after a brief reaction delay. Muscles are then activated, and the driver's lower extremity response will be affected by the muscle active contraction.

According to Behr's experiments [START_REF] Behr | Posture and muscular behaviour in emergency braking: An experimental approach[END_REF], after the driver starts braking, the vehicle will start decelerating at an average acceleration of 0.7 G. Each simulation began while the vehicle hit the barrier and lasted for 120ms. In this paper, we assume that muscle active contraction will be released after the femur force reaches its peak value. (2) In order to compare the femur resultant bending moment in different cases, sixteen cross sections were selected along the femur, referring to the previous study [START_REF] Mo | Incidences of various passenger vehicle front-end designs on pedestrian lower limb injuries[END_REF]. The bending moments of each section at the moment of the peak axial femur force were recorded in Fig. 9.

Results

KTH injury analysis

Most femoral bending moments in active cases are larger than that in inactive cases at the same time and in the same position. The peak bending moments of active models are higher than that of inactive models by about 10-40Nm. As for the 25% offset impacts where the femur fracture the inactive case the maximum value reached 163.6 Nm.

Figure 9. The max bending moments of cross section along the femur

Tibia injury analysis

As shown in Fig. 10, tibia stress distributions of 6 simulation groups were recorded when tibia maximum stress occurred during the impact. Obvious fractures were noted in the lower tibia of both legs in 25% overlap impact, which might be primarily caused by the great invasion of the brake pedal. The fibula fracture at its distal region was also found in C0 case, and both proximal and distal regions of the fibula were injured in C1 case. Besides, the upper right fibula fracture occurred in A1 case. Generally, groups with active muscles show higher stress compared with inactive groups as well as bone injury severity.

Discussion

Effect of active muscles on lower limb injuries was discussed in previous studies with simplified models and loading conditions. A model of the realistic impact environment was established with driver's muscle active contraction, the invasion of the driver's compartment and the response of the restraint system taken into consideration in this study. Driver's muscle active levels during the emergency braking were obtained by volunteer emergency braking experiments.

Considering the drivers' actual response and the repeatability of the driving simulation, there is an assumption in this study defined as follows: when the femur force reaches its peak value, the muscle contraction will be released. The influence of this assumption on the biomechanical responses of lower limbs is non-predicted. The present study used normalized EMG data from experiments for muscle activation level A(𝑡). EMG signals were related to muscle activation but still limited. And in the experiments, we can only measure part of muscles conveniently by surface EMG instrument. The present study only implements muscle activations for four dominant muscles related to braking behaviour. In the future, other methods for defining muscle activations can be studied. Our recent study indicates that combining musculoskeletal model method for pre-defined muscle activation can be a good choice for defining muscle activation in FE human body model [START_REF] Mo | Implementation of controlling strategy in a biomechanical lower limb model with active muscles for coupling multibody dynamics and finite element analysis[END_REF]. In this case, more complicated human movement during driving can be simulated and studied. Additionally, the present study adopts one specific car model and three typical impact conditions. Different kinds of car models and more real impact conditions can be considered for further investigating the robustness of the present findings.

The simulation results show that the reactions of muscles when drivers realize the accident can significantly affect the kinematics and dynamics of lower limbs in the crash. As for KTH injuries, muscle contraction can lead to an obviously additional 1kN axial load on the femur at the moment of the impact. The femur resultant bending moment of active models is also higher than that of inactive models by about 10-40Nm. Similar phenomenon can be found in the previous study [START_REF] Nie | Active muscle response contributes to increased injury risk of lower extremity in occupant-knee airbag interaction[END_REF]. It can be concluded that muscle activity may increase KTH injury risk by increasing both femur axial force and femur bending moment. Such mechanism was also discussed in the research of Chang et al. [START_REF] Chang | Development of a finite element model to study the effects of muscle forces on knee-thigh-hip injuries in frontal crashes[END_REF][START_REF] Chang | Predicting the effects of muscle activation on knee, thigh, and hip injuries in frontal crashes using a finite-element model with muscle forces from subject testing and musculoskeletal modeling[END_REF].

Significantly, fracture of the driver's left femur occurred in both models in 25% offset impacts. The peaks of the femur axial force in models with and without active muscles are 6.8 kN and 5.97 kN, respectively. The values are both lower than 10 kN, which was defined by FMVSS 208 as a femur injury threshold and is now widely used. However, relatively high bending moment of the driver's left femur occurred in 25% offset impact, especially in the active case (175.3 Nm) (Fig. 10). This indicates that the combined action of axial force and bending moment of femur can be a main reason inducing femoral fracture. Thus, the bending moment should be involved in evaluating the femur injury during frontal crashes.

As for tibia injuries, at the initial stage of the impact, the tibial index of the model with muscle activation was about 0.1 higher than that of the model with non-activated muscles (Fig. 11). The active muscle force is considered as a main reason, which may increase both the axial force and the bending moment of tibia [START_REF] Nie | Active muscle response contributes to increased injury risk of lower extremity in occupant-knee airbag interaction[END_REF]. The TI values of the fractured tibias in all cases are far smaller than 1.0, which is widely used as a threshold to evaluate the tibia injury. The fracture occurs even when the TI value is around 0.6. And the tibia fractures primarily occurred around distal end of the tibia in this study. However, the currently widely used TI threshold values were generally based on the experimental results of the tibia mid-shaft. As the previous studies have noted that the tibia tolerance is significantly dependent on the injury area [START_REF] Kerrigan | Experiments for establishing pedestrian-impact lower limb injury criteria[END_REF][START_REF] Mo | Injury tolerance of tibia for the car-pedestrian impact[END_REF][START_REF] Kerrigan | Tolerance of the human leg and thigh in dynamic latero-medial bending[END_REF].

This phenomenon indicates that the TI index can evaluate combined loading but does not take into account the variation of the injury tolerance along the tibia shaft.

Conclusion

This study aims to investigate the influence of active muscles on lower limb injuries with a new-developed lower limb FE model. The frontal impact loading sequence was set with an assumption defined in this paper, which aims to simulate a realistic impact environment.

Besides, driver's emergency braking action, car body structure invasion and response of the restraint system were both taken into consideration. Three types of frontal impacts with different overlap rates were simulated and currently widely used lower limb injury indexes were implemented to evaluate the lower limb injury.

According to the simulation results, the non-ignorable effect of active muscles has been confirmed. The muscle contraction during the impact can significantly affect the kinematics and dynamics of lower limbs by inducing extra loads, which can explain the difference of the injury results between the real accidents and the cadaveric experiments in a way. Additionally, currently widely used lower limb injury criteria can be still lack of robustness. The injury indexes of femur and tibia were both lower than their corresponding thresholds when the injury occurred. Therefore, further research should focus on proposing comprehensive quantified indexes of lower limb injuries considering the active muscle effect, the combined action of different loads and the injury area.

  model geometry was extracted by software Geomagic and Mimics from the computed tomography (CT) scans of bones and Magnetic resonance imaging (MRI) scans of soft tissue. The FE model was developed in LS-DYNA codes by HyperMesh. Considering the influence of element size on computational accuracy and efficiency, most of element sizes were controlled between 4 and 8 mm. As shown in Fig. 1, each muscle of the lower limb-pelvis model consists of three parts: tendon units modelled by non-linear spring elements, muscle passive property modelled by hexagonal elements and muscle active contraction modelled by 1D elements. Muscle passive property was modelled by the Ogden constitutive model in Lsdyna.

Figure 1 .

 1 Figure 1. A lower limb-pelvis FE model with 3D active muscles

Figure 2 .

 2 Figure 2. The frontal impact FE model with driver's emergency reaction for realistic

Figure 3 .Figure 4 .Figure 5 .

 345 Figure 3. Acceleration curves of B pillar

Figure 6 .

 6 Figure 6. Frontal impact loading condition

Figure 7 .

 7 Figure 7. Typical injury occurrences of KTH in 25% overlap impact

Figure 8 .

 8 Figure 8. Femur axial forces in three impacts

Figure 10 .

 10 Figure 10. The peak stresses in the simulation of three impact types (GPa)

Figure 11 .

 11 Figure 11. Tibia index in three impacts

Table 1 .

 1 Muscle active levels A(𝑡)Additionally, bone tissues were modelled with shell and solid elements. According to the anatomy of bone structure, different shell thicknesses and solid densities were used to simulate the epiphysis regions of long bones. The previous study found that such modelling method provides better prediction of strain distribution and bone fracture in the epiphysis regions than the uniform thickness approach[START_REF] Untaroiu | A numerical investigation of mid-femoral injury tolerance in axial compression and bending loading[END_REF]. The whole model consists of 22 muscle groups. Both passive property and active property of this model have been validated against recent relevant experiments. Model validation summary for occupant safety has been shown in the Appendix Table

	Time	Muscle	A(𝑡)%
		musculi tibialis anterior	94±4
		gastrocnemius	102±4
	Start of braking
		rectus femoris	97±5
		biceps flexor cruris	87±4
		musculi tibialis anterior	46±5
	Peak of braking	gastrocnemius	60±5
	force	rectus femoris	36±5
		biceps flexor cruris	26±5

Table 2

 2 Simulation matrix of different frontal impact types and muscle status

		100% full	40% offset	25% offset
		impact	impact	impact
	Inactive muscle	A0	B0	C0
	Active muscle	A1	B1	C1

Table 1 The tolerance of the lower limb model

 1

		Test	Type	Direction	Peak	Tolerance	Reference
		Femur 3-point bending			3.82 kN	3.30.9kN	Yamada, et al.
		Tibia 3-point bending	Quasistatic A-P * : Mid	3.338 kN	2.910.97kN	et al. 1970 [30] Ehler,
		Fibula 3-point bending			0.449 kN	0.442kN	1970 [31] Asang.
							1976 [32]
	Component	Femoral head Static pressure	Quasistatic	Standing posture	7.3 kN	8.43.0kN	Keyak, et al. 1997 [33]
	level				3.945
	validation	Femur 3-point			kN	4.0460.136kN
		bending			398.4	409.84N•m
					N•m		Kerrigan,
					3.15		et al.
		Tibia 3-point bending	Dynamic	L-M * : Mid	kN 233.3	4.041.125kN 25471N•m	2003 [27] Nyquist,
					N•m		et al.
					0.44		1985 [23]
		Fibula 3-point			kN	0.43kN
		bending			31.64	21.7N•m
					N•m
		Foot axial impact	Dynamic	Axial	4.96 kN	5.25kN	Kitagawa, et al. 1998 [34]
	System	Pelvis lateral impact	Dynamic	L-M *	3.55 kN	3.572.03kN	Guillemot, et al. 1997 [35]
	level validation	Knee bolster impact	Dynamic	Axial	9.63 kN	10.10.3kN	Hayashi, et al. 1996 [36]
		KTH knee impact KTH knee impact	Dynamic	Axial	7.13 kN 5.81 kN	7.591.58kN 5.71.38kN	Rupp, et al. 2002 [8]

*A-P: Anterior-Posterior, L-M: Lateral

Acknowledgements

This work is supported by National Natural Science Foundation of China (Grant No. 51621004, 51875187), Hunan Province Science and Technology Plan (Grant No. 2019JJ40021) and Foundation of Zhejiang Province Key Laboratory of Automobile Safety.


Fuhao Mo

Zhefen Zheng

Haotian Zhang

Guibing Li

Zurong Yang

Deyi Sun

In vitro compressive properties of skeletal muscles and inverse finite element analysis: comparison of human versus animals

Keywords: Skeletal muscle, Human gastrocnemius, Inverse analysis, Finite element simulation, Parameter identification

Virtual finite element human body models have been widely used in biomedical engineering, traffic safety injury analysis, etc. Soft tissue modeling like skeletal muscle accounts for a large portion of a human body model establishment, and its modeling method is not enough explored. The present study aims to investigate the compressive properties of skeletal muscles due to different species, loading rates and fiber orientations, in order to obtain available parameters of specific material laws as references for building or improving the human body model concerning both modeling accuracy and computational cost. A series of compressive experiments of skeletal muscles were implemented for human gastrocnemius muscle, bovine and porcine hind leg muscle. To avoid long-time preservation effects, all experimental tests were carried

Manuscript

Click here to view linked References out in 24 hours after that the samples were harvested. Considering computational cost and generally used in the previous human body models, one-order hyperelastic Ogden model and three-term simplified viscoelastic quasi-linear viscoelastic (QLV) were selected for numerical analysis. Inverse finite element analysis was employed to obtain corresponding material parameters. With good fitting records, the simulation results presented available material parameters for human body model establishment, and also indicated significant differences of muscle compressive properties due to species, loading rates and fiber orientations. When considering one-order Ogden law, it is worthy of noting that the inversed material parameters of the porcine muscles are similar to those of the human gastrocnemius regardless of fiber orientations. In conclusion, the obtained material parameters in the present study can be references for global human body and body segment modeling.

Introduction

Finite element (FE) Human body models (HBMs) are widely used for biomechanical analysis in biomedical engineering, transport safety, sport science, and other human involved activities [START_REF] Florio | Effectiveness of various isometric exercises at improving bone strength in cortical regions prone to distal tibial stress fractures[END_REF][START_REF] Halloran | Concurrent musculoskeletal dynamics and finite element analysis predicts altered gait patterns to reduce foot tissue loading[END_REF][START_REF] Hedenstierna | Evaluation of a combination of continuum and truss finite elements in a model of passive and active muscle tissue[END_REF][START_REF] Lenhart | Prediction and Validation of Load-Dependent Behavior of the Tibiofemoral and Patellofemoral Joints During Movement[END_REF]. Skeletal muscle accounts for about 40% of the body mass, therefore an understanding of its mechanical properties is important for virtual human body modeling and related biomechanical application [START_REF] Chomentowski | Skeletal Muscle Mitochondria in Insulin Resistance: Differences in Intermyofibrillar Versus Subsarcolemmal Subpopulations and Relationship to Metabolic Flexibility[END_REF], Takaza, et al., 2013). Muscle tissues can be roughly characterized to be highly nonlinear, viscoelastic, anisotropic and strain rate dependency like most biological tissues, while their properties can also vary significantly due to diversity of species. Thus, the experimental test and numerical modeling of muscle tissues are challenging works.

Skeletal muscle is composed of about 70-80% water, 3% fat and 10% collagen [START_REF] Vignos | A biochemical study of certain skeletal muscle constituents in human progressive muscular dystrophy[END_REF]. As it is very soft, most of the previous studies often tested it in compression setups. Partly due to experimental condition limitation, many previous studies primarily investigated mechanical behaviors of skeletal muscles in compression with different animal species, and obtained corresponding parameters with inverse analysis. For example, several previous studies [START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF][START_REF] Van Loocke | A validated model of passive muscle in compression[END_REF], 2008;[START_REF] Van Sligtenhorst | High strain rate compressive properties of bovine muscle tissue determined using a split Hopkinson bar apparatus[END_REF][START_REF] Song | Dynamic and quasi-static compressive response of porcine muscle [J][END_REF][START_REF] Chawla | Characterization of human passive muscles for impact loads using genetic algorithm and inverse finite element methods[END_REF] investigated high strain rate effects of compression loading on skeletal muscles till 3700/s, while other studies presented experimental results with relatively lower strain rates till 0.001/s [START_REF] Sacks | Biaxial mechanical evaluation of planar biological materials[END_REF][START_REF] Van Sligtenhorst | High strain rate compressive properties of bovine muscle tissue determined using a split Hopkinson bar apparatus[END_REF][START_REF] Song | Dynamic and quasi-static compressive response of porcine muscle [J][END_REF]Xu et al., 2011). With regard to fiber orientations, continuous researches were implemented with different animal muscles primarily considering the fiber orientations of 0°, 45° and 90° (Bosboom et al., 2001 a, b;[START_REF] Van Loocke | A validated model of passive muscle in compression[END_REF][START_REF] Böl | Compressive properties of passive skeletal muscle-The impact of precise sample geometry on parameter identification in inverse finite element analysis[END_REF][START_REF] Mohammadkhah | The in vitro passive elastic response of chicken pectoralis muscle to applied tensile and compressive deformation[END_REF]. Although the uniform macro feature of experimental curves can be noted for different animal muscles, the mechanical properties in detail during compression change a lot due to the species as well as experimental conditions. Thus, the limitations of these animal muscle experiments for FE human body modeling have always existed.

To our knowledge, previous researchers on human skeletal muscle properties under compression loading were not enough and often limited to cadaver specimens with long-time preservation. [START_REF] Vannah | Indentor tests and finite element modeling of bulk muscular tissue in vivo [J][END_REF] conducted quasi-static compression experiments of human muscles to provide design data for lower limb prostheses and other soft tissue supporting devices. [START_REF] Dhaliwal | Structural response of lower leg muscles in compression: a low impact energy study employing volunteers, cadavers and the hybrid III[END_REF] conducted low energy impact tests on human cadaver muscles and obtained corresponding dynamic compressive mechanical properties. These previous studies primarily used the samples harvested from cadavers without indicating details of preservation conditions, like preservation time. However, the previous investigations have indicated that long term formalin preservation can significantly change mechanical properties of biological tissues, even for bone tissues [START_REF] Kikugawa | Effect of long-term formalin preservation on bending properties and fracture toughness of bovine compact bone[END_REF][START_REF] Ohman | The effects of embalming using a 4% formalin solution on the compressive mechanical properties of human cortical bone[END_REF][START_REF] Zhang | The Effect of Formalin Preservation Time and Temperature on the Material Properties of Bovine Femoral Cortical Bone Tissue[END_REF]. Hence, we believe that implementing the experiments with short time preservation samples should be more available for soft tissues like skeletal muscle.

Hence, the objective of the present study was to investigate the compressive properties of human skeletal muscles with short term preservation less than 24 hours and compare them with those of other animals, in order to provide available modeling parameters for FE human modeling or related biomechanical application. In addition, many virtual FE human body or body segment models with skeletal muscle modeling were developed in the previous studies for various considerations, where the Ogden or QLV laws were usually used for these complicated model establishment (Mo et al., 2019). Here, we believe that in a FE human body modeling both accuracy and calculating efficiency should be considered for engineering application. Thus, two generally used material types for the global human body or segment models were selected to obtain specific material parameters by inverse finite element analysis in order to use them in existing or future FE human body models.

Materials and methods

Sample preparation and experimental procedure

All experiments were approved by the ethics committee of the Second Xiangya Hospital of Central South University (NO. 2012-S231) and the informed consent was obtained from the patients participating in the study. Cubic muscle samples with a characteristic edge length of 10 mm were obtained from human gastrocnemius muscle as well as bovine and porcine hind legs. Summary of sample sizes is listed in Table 1.

The bovine samples were harvested from the hindquarter of a large Yorkashire, and the porcine samples were harvested from the hindquarter of a Chinese yellow cattle. The human gastrocnemius muscle was obtained from two amputee donators (age 69±3 years old) without any leg muscle disorders or diseases. The tissues were fixed using formalin to prevent tissue dissolution for no more than 24 hours before the experiments.

Before the experiment, the skin, fat and fascia were firstly removed from the muscle tissues to obtain standard cubic testing specimens. All operations were carried out by a surgeon through a scalpel and an assistive block with ruler to control the size.

Subsequent compressive experiments were performed using different loading rates and along different fiber orientations (0°, 45° and 90°) as shown in Fig. 1. The MTS E45 electronic universal testing machine was used to compress the specimens up to 50% strain at rates of 0.05%/s, 0.5%/s and 5%/s. The measurement range of the selected force sensor is 30kN with a resolution of 0.001N. As pre-compression may lead to unpredictable permanent tissue deformation, each specimen was tested once and no precondition was applied to the specimen. At least three successful specimen tests were achieved at each testing condition. A total of successful 99 compression experiments were conducted at different loading rates and fiber directions.

Mathematical modeling

𝑊 = ∑ 𝜇 𝑖 𝛼 𝑖 𝑛 𝑖=1 (𝜆 1 𝛼 𝑖 + 𝜆 2 𝛼 𝑖 + 𝜆 3 𝛼 𝑖 -3) + 𝐾(𝐽 -1 -𝑙𝑛 𝐽) (1)
Where W is the strain energy of the model; 𝜇 is the initial shear modulus, 𝛼 is the deviatoric exponent, 𝜆 𝑖 is the deviatoric principal stretches in three directions: x, y and z; J represents the bulk modulus. The skeletal muscle was considered incompressible, thus, J equaled to 1. Therefore, the simple one-order Ogden law is expressed as follows,

𝑊 = 2𝜇 𝛼 2 (𝜆 1 𝛼 + 𝜆 2 𝛼 + 𝜆 3 𝛼 -3) (2)
The viscoelastic QLV law is expressed as equation ( 3) and ( 4),

𝜎(𝑡) = ∫ 𝐺(𝑡 -𝜏) 𝑡 0 𝜕𝜎 𝑒 (𝜆) 𝜕𝜆 ⋅ 𝜕𝜆 𝜕𝜏 𝑑𝜏 (3)
Where 𝜕𝜎 𝑒 (𝜆)/𝜕𝜆 is the instantaneous elastic response, 𝜕𝜆/𝜕𝜏 is the stretch history;

𝐺(𝑡 -𝜏) is the relaxation function.

Based on the QLV law, we used a simplified model focusing on instantaneous elastic response in compression. And three terms of model coefficients were adopted. The equation is expressed as follows,

𝜎 𝜀 (𝜀) = 𝐶 1 𝜀 1 + 𝐶 2 𝜀 2 + 𝐶 3 𝜀 3 (4)
Where 𝐶 1 , 𝐶 2 and 𝐶 3 are the coefficients of instantaneous elastic responses; 𝜀 1 , 𝜀 2

and 𝜀 3 are the principal strains in three directions.

Inverse analysis

Inverse analysis has been applied for defining material parameters of soft tissues ( 2005,2007). The purpose of the inverse analysis was to minimize the error between the experimental results and simulation responses by using an optimization algorithm to identify optimal material parameters [START_REF] Böl | Compressive properties of passive skeletal muscle-The impact of precise sample geometry on parameter identification in inverse finite element analysis[END_REF]. The agreement between experimental results and simulation results were evaluated by the objective function as follows,

n n i F i F i x f FEM    ) ( exp 2 ) ( (5)
Where 𝐹 𝑖 𝑒𝑥𝑝 is the experimental force value, 𝐹 𝑖 𝐹𝐸𝑀 is the responsive force of finite element simulation, and n is the number of data points. The process of inverse analysis for material parameter identification is described in the following section.

Identification of material parameters

To identify material parameters of the tested skeletal muscles, a method combing inverse solve and finite element optimization of material parameters were employed, as shown in Fig. 2. The errors between the simulation and experimental results are calculated by equation [START_REF] Chawla | Characterization of human passive muscles for impact loads using genetic algorithm and inverse finite element methods[END_REF]. A custom python program using a least-squares algorithm was firstly used to fit experimental results to acquire initial parameters of the corresponding laws. Then, their values were optimized and finally determined by inversed finite element analysis. Adaptive response surface method (ARSM) was adopted in the FE simulations to minimize the error and obtain optimal parameters until calculation convergence referring to the previous study (Wang et al.,2007), which is a direct optimization method based on response surface. The basic idea of ARSM method is to construct a response surface with fewer sample points first and determine the optimization direction. During the optimization process, new design points are obtained along the gradient direction of the response surface, and introduced into the design space so that the response surface model can be gradually updated to improve optimization efficiency and fitting accuracy. By optimizing and updating the response surface, the optimal solution can be continuously approached to the real optimal solution.

Fig.2. Schematic illustration of material parameter identification

The FE simulations were implemented with Ls-Dyna codes. The specimen was meshed with hexahedral solid elements of the size 1 mm by Hypermesh software. Element quality was checked to avoid distorted elements and to satisfy the minimum time step criteria. The time step was set as 1.0e-5. Surface to surface contacts were defined between the platens and the specimen, and a self-contact of the specimen was defined to avoid solid collapse during large deformation. The lower platen was constrained with no freedom in all directions, while the upper platen can freely translate along the compression direction. The friction coefficient between the platen and specimen was set as 0.3 according to the previous study [START_REF] Wu | Analysis of effects of friction on the deformation behavior of soft tissues in unconfined compression tests[END_REF]).

Results

Experimental results

Characteristic compressive stress-strain responses are classified based on tissue fiber orientations, compression strain rates and species (Fig. 3-Fig. 5). The bovine and porcine muscle properties are illustrated in Fig. 3 shown in Fig. 6, the substantial difference can be noted due to species, fiber orientations and strain rates. The compressive properties of the porcine muscles seem to be closer to the human gastrocnemius muscles than those of the bovine muscles.

In most testing conditions except for 90° bovine muscle, potential strain effects are noted in the testing range from 0.05% to 5% for all species and fiber orientations, while the Cauchy stress value inclines significantly with the increasing strain rate at 50% strain as shown in Fig. 6 (p<0.005). With the increase of the strain rate, the scattering range of the compressive curves in most of testing conditions tends to enlarge. The stiffest tissue response was obtained when the muscle fibers were oriented in perpendicular to the loading direction (90° fiber orientation), followed by a fiber orientation of 45°. Fiber orientation of 0° in parallel with the loading direction yielded the softest response.

Identification of material parameters

We used the averaged experimental curves of different fiber orientations and different species as the objective curves for the optimization of material parameters. Through the above-mentioned method, the typical fitting results of numerical simulations to experimental results considering the human skeletal muscles of three fiber orientations are shown in Fig. 7. All simulation results are in good agreement with the experimental responses.

Inversed material parameters of different muscle tissues using one-order Ogden law were illustrated in Table 2. The corresponding material parameters of three-terms simplified QLV law were listed in Table 3. For both material laws, the obtained material parameters for different species and fiber orientations present some differences.

Although significant differences (p< 0.05) of material parameters existed between the porcine and human muscles considering different fiber orientations, the inversed material parameters of the averaged curves of all three fiber orientations show similar values for both human (μ=3.43 kPa, α=8.74) and porcine (μ=3.63 kPa, α=8.74) muscles concerning one-order Ogden law.

Discussion

Although The Cauchy stress-strain relationship presented in this study provided a reference for investigating mechanical behaviors of human skeletal muscles. The results showed that the Cauchy stress during compression increased with the increasing strain rate, which was similar to the previous researches (Chawla, et 2008). The cross-fiber orientation mechanical response was found to be stiffer than other fiber orientations, which was similar to the findings of the previous researches with different loading strain rates [START_REF] Song | Dynamic and quasi-static compressive response of porcine muscle [J][END_REF][START_REF] Böl | Compressive properties of passive skeletal muscle-The impact of precise sample geometry on parameter identification in inverse finite element analysis[END_REF].

As the inversed material parameters shown in Table 2, the optimized shear modulus μ The present study selected one-order Ogden law for skeletal muscle modeling, as it is a simple and reasonable solution for soft tissue modeling. Especially concerning its application in a global biomechanical model, the high computational efficiency was extremely appreciated in our previous foot model and lower limb modeling [START_REF] Mo | A Lower Limb-Pelvis Finite Element Model with 3D Active Muscles[END_REF](Mo et al., , 2019)). We compared the computational cost of one-order Ogden law and threeterms simplified QLV law in the present compression modeling using Ls-Dyna solver.

With 8-core Xeon® CPU, the calculation time for a simple compression simulation is 23 min and 27mins, respectively. The one-order Ogden law has some advantages in saving computational cost as well as controlling the stability of a global human body model for loading conditions without large concerns on viscoelasticity. Considering a global human body or large body segment modeling, we believe the balance of modeling accuracy and computational cost should be achieved properly with regard to research objectives.

The present study combined experimental tests and inverse analysis method to investigate the effects of loading rate, fiber orientation and species differences on compressive behaviors of skeletal muscles. A large number of experiments and numerical simulations were implemented, and available material parameters were obtained for human body modeling reference. There are still several limitations. First, the limited sample size for each testing condition may partly limit part of the robustness of our results. As limited pure muscle tissues can be obtained after removing other fiber tissues in one donator, we had to harvest samples from two amputees not as we can obtain all muscle samples from one animal. But if we enlarge our sample size, the effects of individuals' difference would be also extended. Second, we use an ideal and standard geometry for inverse finite element analysis to obtain muscle material parameters. As indicated by the previous study [START_REF] Böl | Compressive properties of passive skeletal muscle-The impact of precise sample geometry on parameter identification in inverse finite element analysis[END_REF], to further precisely define the muscle parameters, an inverse finite element method with detailed geometry definition could be more powerful. Additionally, we also believe that using short-time preservation tissues for experiments or adopting in-vivo measuring methods like the SWE method for obtaining or evaluating material parameters can be further investigated for virtual FE human model establishment.

Conclusion

Inspired by the objective of human body modeling, the present study investigated the compressive properties of skeletal muscles with different species, loading rates and fiber orientations. Two commonly used material laws of Ogden and QLV constitutive models were selected and simplified for numerical analysis in the consideration of computational efficiency and modeling accuracy. Optimized material parameters were well obtained through inverse analysis based on custom-made codes and optimal finite element simulations, which can be available references for human body modeling. The experimental results show substantial differences between species, loading rates and fiber orientations. Referring to the material parameters of one-order Ogden law, similar values between the porcine muscles and the human gastrocnemius were found regardless of fiber orientations. The inversed shear modulus was further verified by in vivo SWE measurement, thus it also indicated that the present parameters can reflect the properties of in vivo human gastrocnemius muscle. Coupling experimental and simulation investigations on different human muscles could be extended in the future to further evaluate the present method and findings.

Fig. 1 .

 1 Fig.1. Schematic of compression experiments with different loading rates and fiber directions

Following

  the previous studies (Comley et al. 2012, Takaza et al. 2013, Mo et al, 2018, 2019a, b), one-order Ogden hyperelastic material law and three-terms quasi-linear viscoelastic (QLV) material law were usually selected for modeling passive behavior of muscle tissues in human body models. The Ogden hyperelastic material law is expressed as follows,

  and Fig. 4, respectively. Compressive responses of the human gastrocnemius muscles are shown in Fig. 5. As general soft biological tissues, the experimental data of the muscle tissues scatter on a large scale in the present study. When comparing compressive responses of different species as

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 3456 Fig.3. Compressive stress-stretch responses of the bovine muscles at various rates and fiber orientations

Fig. 7 .

 7 Fig.7. Typical fitting results of the finite element simulations to the experimental results of the human skeletal muscles.

  some previous studies have investigated skeletal muscle properties in different loading types using different animal specimens (Mohammadkhah et al.,2016, Nie et al.,2011, Simms et al.,2012, Takaza et al.,2013a), the experimental results of human skeletal muscles especially the data with short preservation time in formalin or in status were rare. In the present study, we aim to provide material properties for human body modeling, especially in the application of comfort and impact analysis during seating and car crash. In these two loading conditions, human body muscles generally experience a compression loading. In addition, when we compared the published muscle properties of different animals (Song et al., 2007, Dhaliwal et al., 200, Mohammadkhah et al.,2016, Van Sligtenhorst et al.,2006, Takaza et al.,2014), substantial differences of mechanical behaviors in compression can be noted. Hence, the first purpose of this study was to investigate compressive behaviors of human skeletal muscle without long formalin storing time, then obtain their material parameters of generally used and computationally efficient material laws for the improvement or establishment of human body models. By the way, the dominant fiber directions can be different considering anatomic structures of different muscles. Muscle compressive properties from three representative fiber directions were measured and analyzed. Compressive properties of muscle tissues in other fiber directions should be located between them. In the establishment of a related human body model, we recommended either establishing or employing an anisotropic constitutive material law in the model, or optimizing the corresponding muscle parameters of the available isotropic laws considering fiber distribution and boundary loading conditions.

of 5 .

 5 07 kPa for human gastrocnemius was close to in vivo experimental results using supersonic shear wave elastography (SWE) of Maïsetti et al. (2012) study or the results of our verifying measurements. This can indicate the reliability of our results, and that the obtained tissue properties are close to in vivo muscle tissue properties. In addition, the average inversed material parameters of one-order Ogden law regardless of fiber orientations for the porcine and human muscles are similar. But those of the bovine muscle are largely different with them. To some extent, it could indicate that using porcine muscle for experimental tests could provide more useful information for human muscle tissue modeling if the human tissue experiment is limited. But it can be also influenced by native species' difference. To extend experimental possibility, we believe that it is important to find proper alternative animal tissues with similar mechanical properties to human tissues for experimental researches.

Figure 1 .

 1 Figure 1. Schematic of compression experiments with different loading rates and fiber directions.

Figure 2 .

 2 Figure 2. Schematic illustration of material parameter identification.

Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 3456 Figure 3. Compressive stress-stretch responses of the bovine muscles at various rates and fiber orientations. (a) 0° (b) 45° (c) 90°

Figure 7 .

 7 Figure 7. Typical fitting results of the finite element simulations to the experimental results of the human skeletal muscles. (a) 0° (b) 45° (c) 90°

Figure. 3

 3 Figure.3(a)

Table 1 . Sample size and average specimen dimensions at different fiber orientations

 1

	Specimen	Fiber direction	Average height ±std(mm)	Average width ±std(mm)	Average thickness ±std(mm)	Sample size (N)
		0°	10.5±0.6	10.7±0.5	10.5±0.6	12
	Bovine	45°	10.8±1.2	11.1±1.2	10.6±0.6	12
		90°	10.9±0.9	10.8±0.7	10.8±0.7	12
		0°	11.3±0.9	10.8±0.8	10.9±0.9	12
	Porcine	45°	11.1±1.4	10.5±1.1	10.5±0.5	12
		90°	11.7±0.8	11.3±0.9	10.8±0.7	12
		0°	10.0±1.0	10.7±0.5	10.3±0.5	9
	Human	45°	9.9±0.9	10.4±0.8	10.1±0.7	9
		90°	12.0±2.1	10.6±1.9	10.3±0.9	9

Table 2 Optimized muscle material parameters of one-order Ogden law

 2

	Species	Fiber orientation	Material parameters	R 2 (FE)	RMSE
		0°	μ=0.26 kPa α=14.00	0.994	1.21
		45°	μ=0.65kPa α=13.83	0.998	1.67
	Bovine
		90°	μ=2.00 kPa α=10.11	0.997	8.91
		Average	μ=1.18 kPa α=11.97	0.998	2.12
		0°	μ=0.54 kPa α=14.00	0.983	11.07
		45°	μ=3.80 kPa α=9.00	0.969	115.50
	Porcine
		90°	μ=3.96 kPa α=8.92	0.994	23.08
		Average	μ=3.63 kPa α=8.74	0.983	40.59
		0°	μ=2.24 kPa α=8.97	0.989	15.72
		45°	μ=3.77 kPa α=9.00	0.982	63.88
	Human
		90°	μ=5.07 kPa α=9.14	0.994	43.79
		Average	μ=3.43 kPa α=8.74	0.934	173.33

Table 3 Optimized muscle material parameters of three-terms QLV law

 3

	Species	Fiber orientation	Material parameters	R 2 (FE) RMSE
		0°	C1=88.61 C2=-772.05 C3=1716.00	0.937	18.74
		45°	C1=169.70 C2=-1500.86 C3=3607.17	0.957	45.73
	Bovine
		90°	C1= 132.00 C2= -1140.00 C3=3320.00	0.960	79.09
		Average	C1= 122.18 C2= -1065.46 C3=2808..00	0.978	26.77
		0°	C1=90.00 C2= -1088.22 C3=2849.00	0.977	22.26
		45°	C1=99.00 C2= -994.16 C3=3366.99	0.992	22.51
	Porcine
		90°	C1=189.00 C2= -1477.87 C3=4488.00	0.991	36.25
		Average	C1= 132.00 C2=-990.00 C3=2981.00	0.982	40.32
	Human	0°	C1= 115.36 C2= -965.06 C3=2758.01	0.967	42.45

  [START_REF] Chawla | Characterization of human passive muscles for impact loads using genetic algorithm and inverse finite element methods[END_REF][START_REF] Van Sligtenhorst | High strain rate compressive properties of bovine muscle tissue determined using a split Hopkinson bar apparatus[END_REF][START_REF] Song | Dynamic and quasi-static compressive response of porcine muscle [J][END_REF][START_REF] Van Loocke | A validated model of passive muscle in compression[END_REF][START_REF] Van Loocke | Viscoelastic properties of passive skeletal muscle in compression: Stress-relaxation behaviour and constitutive modelling[END_REF]. The Cauchy stress of 50% compression strain at the compression loading rate of 5%/s was about 3 times larger than that at the compression rate of 0.05%/s, which was qualitatively similar to the results at other strain rates[START_REF] Van Loocke | A validated model of passive muscle in compression[END_REF] Van Loocke, et al.,

Acknowledgments

This study was supported by the National Key Research and Development Program of China (2018YFC1105800), National Natural Science Foundation of China (Grant No. 51875187, 51621004), and Hunan Province Science and Technology Plan (Grant No. 2019JJ40021). We also thank gratefully to Dr. Tang Liu's and Dr. Pengcheng Dou's help for experimental tests.

Conflict of interest statement

The authors assert that there are no conflicts of interest of any type.

Conflict of Interest Statement

Fuhao MO, Zhefen ZHENG, Haotian ZHANG, Guibing LI, Zurong YANG, Deyi Sun, J Biomechanics Title: In vitro compressive properties of skeletal muscles and inverse finite element analysis: comparison of human versus animals The authors assert that there are no conflicts of interest of any type. This assertion is also included in the manuscript.

Conflict of Interest Statement