Huy Do

Pascal Bourdon

David Helbert

Mathieu Naudin

Remy Guillevin

7T MRI super-resolution with Generative Adversarial Network

The high-resolution magnetic resonance image (MRI) provides detailed anatomical information critical for clinical application diagnosis. However, high-resolution MRI typically comes at the cost of long scan time, small spatial coverage, and low signal-to-noise ratio. The benefits of the convolutional neural network (CNN) can be applied to solve the super-resolution task to recover high-resolution generic images from low-resolution inputs. Additionally, recent studies have shown the potential to use the generative advertising network (GAN) to generate high-quality super-resolution MRIs using learned image priors. Moreover, existing approaches require paired MRI images as training data, which is difficult to obtain with existing datasets when the alignment between high and low-resolution images has to be implemented manually.

This paper implements two different GAN-based models to handle the super-resolution: Enhanced superresolution GAN (ESRGAN) and CycleGAN. Different from the generic model, the architecture of CycleGAN is modified to solve the super-resolution on unpaired MRI data, and the ESRGAN is implemented as a reference to compare GAN-based methods performance. The results of GAN-based models provide generated highresolution images with rich textures compared to the ground-truth. Moreover, results from experiments are performed on both 3T and 7T MRI images in recovering different scales of resolution.

Introduction

In recent years, there has been a huge increase of computer vision research leading to an important proportion of information related to this field. Image superresolution (SR) is an exploitable image processing field that refers to recovering high-resolution (HR) images with rich details from given low-resolution (LR) images.

At the present time, several deep learning-based methods have been proposed to improve the quality of images. Convolutional neural networks (CNN) can model various complex structures and utilize large quantities of training data. After completing the training phase, no additional process is required to produce output images; thus, it does not cost time for the sub-sequence process. Recently, many CNN-based methods have out-performed in the field of SR. Dong et at. [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF] proposed the SRCNN, which learns an end-to-end mapping from LR image to HR images. Since the work of Dong et al. [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF], different CNN based super-resolution approaches have been implemented.

By improving the network architectures and training process, CNN-based SR has achieved significant success on both objective (peak signal-to-noise ratio -PSNR) [START_REF] Lim | Enhanced deep residual networks for single image super-resolution[END_REF][START_REF] Ledig | Photo-realistic single image super-resolution using a generative adversarial network[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF] and subjective (human visual quality assessment tests) [START_REF] Ledig | Photo-realistic single image super-resolution using a generative adversarial network[END_REF][START_REF] Mehdi | Enhancenet: Single image superresolution through automated texture synthesis[END_REF] criteria. A deep neural network with large capacity can improve the performance in general [START_REF] Sun | On the depth of deep neural networks: A theoretical view[END_REF]. However, it also increases the computation cost of the training phase. With different techniques such as skip connection [START_REF] He | Deep residual learning for image recognition[END_REF], embedding [START_REF] Szegedy | Inception-v4, inception-resnet and the impact of residual connections on learning[END_REF] or normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF], reducing computation time for a effective training now is possible. Following the study of Kim et al. [START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF], a deep network using advanced techniques to increase network depth could ease the difficulty of training networks, and the structure of the neural network is the key to obtain the high quality of SR outputs.

Tong et al. [START_REF] Tong | Image super-resolution using dense skip connections[END_REF] have proposed the densely connected network (SRDenseNet) to solve SR tasks, using several hierarchical features in different layers into the final construction layer. Sub-sequence layers can effectively utilize extracted feature information from convolutional layers. It is a new measure to preserve the features of different levels extracted from different convolutional layers in the network. It retains the ground truth image features to a greater extent, so the dense connection can effectively improve image reconstruction quality. SRDenseNet has significantly improved performance over the model using multi-level features, indicating that level fusion is indeed beneficial for SR problems. Besides, since ResNet has been published by He et al. [START_REF] He | Deep residual learning for image recognition[END_REF], it has proved its efficiency in reducing training computation time and benefit in SR fields [START_REF] Gross | Training and investigating residual nets[END_REF]. Since their introduction in [START_REF] Goodfellow | Generative adversarial nets[END_REF], generative adversarial networks (GAN) have become a trend to solve many different computer vision problems. The generic GAN model includes a generator and a discriminator. While the generator is trained to learn a mapping from source images in a domain to target images in another domain, the discriminator distinguishes and targets generated images with binary labels. Once well trained, GAN can perform on high-dimensional target images.

Following the progress of neural networks applied in natural images, medical image analysis is one of several applications that is benefited from the improvement of SR approaches. Further research has been proposed to improve medical image quality, such as computed tomography (CT) or Magnetic Resonance Imaging (MRI).

MRI is widely used in medical imaging because of its non-invasive assessment of the body's anatomy and physiology in both health and disease while providing the best contrast resolution on soft tissues. However, MRI images are normally acquired with a finite resolution limited by the signal-to-noise ratio (SNR), hardware, or time limitation. In the clinical and research centers, usually high-resolution and high-contrast MRI is preferred because it provides critical structural details with a smaller voxel size. Therefore, the request for image quality with sufficient details is fast increasing in medical imaging. In recent years, the ultra-highfield MR 7T scanners were introduced, which allow high-resolution MRI scanning with many advantages compared to routine MRI, such as 3T MRI or 1.5T. The 7T MRI scanner can provide higher resolution images while maintaining a similar signal-to-noise ratio (SNR), which has a linear relation to the magnetic field in general for MRI scanning.

Recently, many studies have been proposed for medical image analysis, especially super-resolution. the variety of methods stretches from statistical method such as interpolation [START_REF] Thomas Martin Lehmann | Survey: Interpolation methods in medical image processing[END_REF][START_REF] Thévenaz | Interpolation revisited [medical images application[END_REF], dictionary mapping [15, 16], self-learning [START_REF] José V Manjón | MRI superresolution using self-similarity and image priors[END_REF] to automatically techniques such as CNN based-method [START_REF] Chen | Brain MRI super resolution using 3d deep densely connected neural networks[END_REF][START_REF] Zeng | A comparative study of CNN-based superresolution methods in MRI reconstruction and its beyond[END_REF] or hybrid methods [START_REF] Qu | Synthesized 7T MRI from 3T MRI via deep learning in spatial and wavelet domains[END_REF]. However, most existing approaches are implemented on private datasets because algorithms require paired datasets for training. Such datasets usually come from collaborating hospitals. However, these are not large enough due to the tedious and time-consuming task of paired data retrieval. As an example, to obtain paired 3T and 7T MRI, patients must take the test on the same day, and later, images have to be aligned manually. Besides, there is also the difficult to publish due to security issues.

For super-resolution task, a paired dataset contains both low and high-resolution MRI, which is not easy to obtain. There is no available scanner to produce HR MRI, which generally have HD or 2K resolution. All the pre-process of super-resolution has to be prepared manually with a very complicated process of alignment. At this moment, there is no public dataset with paired MRI.

This paper implements two different types of GANbased models to solve SR on MRI images: the Enhanced Super-Resolution GAN (ESRGAN) and Cycle-GAN. The advantage of the approach is that image pairing is not required for the training dataset. Therefore they can be executed on several publicly available MRI datasets, thus overcoming the limitations explained earlier.

CycleGAN is known as a very popular method for image-to-image translation. However, the superresolution problem requires some special characteristics due to the difference between input and output. With some modifications to network architecture, we make CycleGAN fit with the SR problem. The generator is modified with several different building blocks from the generic CycleGAN to automatically solve up-sample and down-sample images by applying skip connection techniques. In term of ESRGAN [START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF], it is an improved version of SRGAN [START_REF] Ledig | Photo-realistic single image super-resolution using a generative adversarial network[END_REF]. With these special modifications of network architecture, ESRGAN has shown excellent performance in SR fields. To our knowledge there is no ESRGAN implementation for medical images at the present time, therefore we decide to explore ESRGAN performance in MRI super-resolution. It also can be used as a reference to compare with CycleGAN performance later.

In the next section, we introduce background knowledge and concepts of image super-resolution with different CNN-based methods and techniques. Then we discuss the details of network architectures used in our experiments.

Related work 2.1 Image super-resolution

Before deep learning-based approaches achieved stateof-the-art performance, SR techniques mostly relied on statistical analysis, interpolation, edge-preservation, and sparse dictionary learning [START_REF] Yang | Single-image super-resolution: A benchmark[END_REF]. Since the study of Dong et al. [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF], CNNs have been increasingly popular to solve the SR problem. The SRCNN can handle feature extraction, feature space building and image reconstruction together in an end-to-end training. Later, many follow-up approaches with the improvement of network structures have been inspired by the SRCNN [START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF][START_REF] Tong | Image super-resolution using dense skip connections[END_REF][START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF].

Residual, dense blocks and skip connection

The performance of deep learning methods is fast improving. In contrast, these models are significantly increasing model size, wherein the depth of modes becomes a practical problem. Moreover, when the network is too deep, gradient disappearance and gradient explosion issues are declared [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. Although these gradient issues can be solved by data regulation and batch normalization, it can lead to model performance degradation. To address this problem, ResNet proposed by He et al. [START_REF] He | Deep residual learning for image recognition[END_REF] introduced residual learning, where the output of the previous convolutional layer connected to the next for smoother information flows through a short cut. The detail of residual blocks is shown in Figure 1. The shortcut connection neither increases the number of network parameters nor the computational complexity of algorithms. DenseNet [START_REF] Huang | Densely connected convolutional networks[END_REF] proposed a connectivity pattern to improve the flow of feature information by concatenating information of previous layers. The network is more efficient and outperforms ResNet with fewer parameters. In addition, recent studies in image super-resolution show that removing BN layers in residual blocks [START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF] and dense blocks [START_REF] Lim | Enhanced deep residual networks for single image super-resolution[END_REF] can reduce computational cost, memory usage, and boost model performance. The flow of the gradient is unobstructed due to the direct link between layers. SRDenseNet [START_REF] Tong | Image super-resolution using dense skip connections[END_REF] is a version of DenseNet to solve the SR problems.

Additionally, Zhang et al. [START_REF] Zhang | Residual dense network for image super-resolution[END_REF] introduced Residual Dense Network (RDN), a combination of DenseNet and ResNet to solve the image super-resolution task by using residual dense blocks (RDB), which contain several dense connected layers, a local feature fusion, and the residual learning as the final layer. Figure 1 illustrates the design of all blocks in different superresolution models with the removing of BN layers .

Generative adversarial network

and GAN-based Super-resolution Most previous super-resolution methods aim to optimize HR image reconstruction by minimizing the voxel-wise difference between original and generated images. However, a drawback that merely cares about local pixel-wise differences leads to extreme difficulty in restoring important small details. In contrast, if global perceptual constraints can be taken into account, the SR model is guided by both local intensity information and patch-wise perceptual information, probably leading to a better and sharper SR reconstruction [START_REF] Ledig | Photo-realistic single image super-resolution using a generative adversarial network[END_REF].

With benefits from the GAN framework of Goodfellow et al. [START_REF] Goodfellow | Generative adversarial nets[END_REF] for its unsupervised-learning potential of capturing perceptually important image features, Ledig et al. [START_REF] Ledig | Photo-realistic single image super-resolution using a generative adversarial network[END_REF] proposed the SRGAN to handle the superresolution issue. The principle of GAN relies on two components: a generator and a discriminator. While the generator generates fake images as real as possible, the discriminator distinguishes and evaluates to improve generated image accuracy. In the end, the discriminator can separate real and generated image, while the generator can produce realistic-looking images. Next, the perceptual loss is defined in the GAN model [START_REF] Goodfellow | Generative adversarial nets[END_REF] and extended in SRGAN. SRGAN defines the perceptual loss on the activation layers of a pretrained deep network, where the distance between two activated features is minimized.

Several techniques have been implemented to provide different architectures of building units to transform LR into HR image and reduce computation cost during the training phase. Residual blocks (RB) from ResNet of He et al. [START_REF] He | Deep residual learning for image recognition[END_REF] and Dense block (DB) from DenseNet of Huang et al. [START_REF] Huang | Densely connected convolutional networks[END_REF] are the most popular architecture use in SR tasks. Based on different models, these building blocks can be combined or modified. In SRGAN, residual blocks have been used for the generator.

ESRGAN

Improved from SRGAN, ESRGAN proposed by Wang et al. [START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF] aims to increase the quality of HR images. While SRGAN uses residual blocks as basic blocks, ES-RGAN uses Residual-in-Residual Dense blocks (RRDB) without Batch Normalization (BN) layers as building units. As shown in Figure 1, RRDB is a combination of multi-level residual blocks and dense block connection.

Recent studies have shown that removing BN layers can reduce computational cost, memory usage, and boost model performance [START_REF] Lim | Enhanced deep residual networks for single image super-resolution[END_REF]. Although BN layers use mean and variance computations to normalize the features during training and later during testing, BN layers tend to produce unpleasant artifacts and limit generalization ability when the difference between training and testing set is significant [START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF].

The addition of layer and connections can improve model performance [START_REF] Zhang | Residual dense network for image super-resolution[END_REF][START_REF] Zhang | Image superresolution using very deep residual channel attention networks[END_REF], RRDB exploits deeper and more complex connections than the residual blocks of SRGAN. Meanwhile, the general high-architecture of the ESRGAN model is kept as SRGAN. The details of the generators with the additions of RRDB are shown in Figure 2.

The discriminator of ESRGAN is also improved from SRGAN, based on the relativistic GAN [START_REF] Jolicoeur-Martineau | The relativistic discriminator: a key element missing from standard gan[END_REF]. Instead of estimating the probability of an input being real or natural, the relativistic discriminator calculates the probability of a real image to be relatively more realistic than a fake image [START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF].

In terms of perceptual loss, ESRGAN uses features before the activation layers, which helps overcoming the drawbacks of the original design. The perceptual loss in SRGAN can cause inconsistent reconstructed brightness compared with the ground-truth image or the sparsely activate of features when deep network is very deep [START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF]. The architecture of ESRGAN is implemented from the original as a reference to compare performance of GAN-based model.

CycleGAN

The core problem of the SR task in the real world is the lack of paired datasets. Zhu et al. [START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF] introduced Cycle-GAN -an image-to-image translation framework using unpaired data. Later, it has inspired many following studies in different computer vision tasks, and potential for SR. CycleGAN aims to translate input images from a class into another class without the requirement of Like the RRDB, RDB also removes BN layers to reduce computational time, memory usage, and to speed up the training process. With the addition of the basic block, the network also introduces contiguous memory [START_REF] Zhang | Residual dense network for image super-resolution[END_REF], which allows direct access of information between each block layer. The model also improves the preservation of information inside the network using local and global feature learning which also combine the low and high-level features.

Using the RDN, we focus on the generator architecture of CycleGAN to adapt to super-resolution problems. In the super-resolution CycleGAN, the generator contains several RDB for feature extractions with local and global residual learning to keep information convection. Next, local and global features are stacked for up-sampling or down-sampling to produce output. The details of the generators with the additions of RDB is shown on Figure 3.

A generic CycleGAN model can perform translations between two classes in both directions. Therefore we build two different generators to up-sample from LR to HR and down-sample from HR to LR simultaneously. The architecture of the discriminator is kept to evaluate the performance of generators.

The cycle architecture or CycleGAN is kept identical to the original model [START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF]. The generator should eventually be able to trick the discriminator about the authenticity of its generated images. This can be done if the recommendation factor of the discriminator for generated images is as close to 1 as possible. The content of cycle consistency and adversarial loss -the outstanding properties of CycleGAN -is kept. The adversarial loss of CycleGAN is calculated based on the MSE loss between the real and fake images of two classes. To handle the generating of non-related output from input, the cycle loss consistency is applied. In detail, when an image is generated, the second generator will convert it back to the original class. Then, the cycle consistency is calculated by the average of L1 loss between the real input and the cyclic output. The difference between the original image and the cyclic image should be as small as possible.

The discriminator must be trained such that recommendation factor for images from the first category be as close to 1 as possible, and vice versa for the second discriminator. So the first discriminator would like to minimize the loss between the generated and the real ones, and the same goes for the second category as well. Since discrimination should be able to distinguish between generated and original images, it should also be predicting 0 for images produced by the generator.

Experiment and Results

Dataset

At this moment, there is no benchmark MRI dataset for the super-resolution task. Current MRI scanners are not available to produce high-resolution MRI, which generally have HD or 2K resolution. Obtaining highresolution MRIs requires them to be processed manually through a complex alignment system.

As mentioned in the previous section, we aim to implement methods that does not require paired images as training data. In this project, we used the BraTS2018 [START_REF] Bjoern H Menze | The multimodal brain tumor image segmentation benchmark (brats)[END_REF] dataset -a public dataset containing 3T MRI images with different types of sequences. BraTS is a very popular dataset used in the segmentation and classification domain. The variety of the dataset is confirmed when samples were acquired with different clinical protocols and various scanners from multiple institutions. The field of view is 155 × 240 × 240, wherein the slice thickness is 1mm. The size of the original image is kept at a 1:1 ratio when the output is 240 × 240 pixels.

Within the project scope, we also want to explore the performance of methods on different types of MRI. 7T MRI has proved its function in medical image analysis. However, there is no 7T dataset available at this moment. With the support of the 7T system from Siemens Healthineers at Poitiers University Hospital, we can build a dataset containing several brain 3T and 7T MRI samples with high-quality. We aim to use that dataset to explore super-resolution performance. Different from 3T, 7T technology provides MRI images with higher quality and deeper voxel spacing between slices. The field of view is 255 × 320 × 320, wherein the slice thickness is 0.7mm. The original image's size is also kept at a 1:1 ratio when the output is 320 × 320 pixels.

Training setting

The training process is implemented on BraTS2018, which contains only 3T MRI. Later, the testing set of BraTS and 7T MRI dataset is used on the testing phase to evaluate the performance of the superresolution model. Due to the lack of high-resolution (HD or 2K) MRI, training images in the dataset are degraded into lower-resolution to simulate LR images with scale factors. All image resolution changes are implemented during the pre-processing step without modifying original data during the training phase.

The training of both ESRGAN and CycleGAN is executed on patches to ensure the diversity of data. Besides, we also want to reduce the computational cost of the model during the training process. The complexity of GAN-based models is very considerable, with millions of parameters. This number will increase along with the increase of model depth or the size of the input. For each training batch, images are randomly extracted into 16 patches on LR images corresponding to a size of 64 × 64 on HR images.

For CycleGAN, we keep the content of adversarial loss and cycle consistency loss. The generator is modified to contain residual dense blocks. In CycleGAN, due to the network architecture, two generators are built to up-sample and down-sample images. In the downsample generator, there are two options to down-sample images with convolutional layers or pooling layers. The learning rate is initialized to 0.0005 and decay starts after every 100 iterations. The ADAM optimizer is used to update network weights based on training data.

For ESRGAN, we have implemented the model following the original paper to use it as a reference for the GAN-based method reference. The generator contains several RRDB. Different from CycleGAN, a single ESRGAN model only can implement the upsample or down-sample images. The learning rate is initialized to 0.0002 and decay starts after the fifth epoch. The ADAM optimizer is also used to update network weights. The first 100 batches are to "warm-up" the training phase, which uses only pixel-wise loss.

Results

Benefiting from these improvements, results of GANbased models provide generated high-resolution images with rich textures that compare to ground-truth data. Results from our experiments are performed on both 3T and 7T MRI images to reconstruct different resolution levels. Figure 4 1 shows average distortion and fidelity values obtained on the test set by GAN-based methods. We have also included simple bicubic interpolation to use it as a reference for quality. Image quality measurement in terms of PSNR and SSIM shows that CycleGAN methods can achieve relatively low distortion at the ×4 scale factor, while ESRGAN is still limited. We consider that the CycleGAN has better and more stable performance compared to ESRGAN.

As shown on Figure 4 and Table 1, both ESRGAN and CycleGAN provide super-resolution 3T MRI images with rich textures. The performance of CycleGAN outperforms that of ESRGAN especially on 7T MRI data. CycleGAN SSIM scores exceed 0.9, thus ensuring a better fidelity of image structures such as contours or fine details compared to the original images. The generated images of ESRGAN contains noise, while the CycleGAN images are relatively stable. We also test the performance of models on 7T MRI data.

The CycleGAN model is a prevalent method known for image synthesis. With the advantage of self-learning between two different classes, it can be used for superresolution tasks. By the modification of network architecture, we use CycleGAN to match the problem of SR.

To have a reference to compare with CycleGAN, we also implemented ESRGAN -a popular method in the super-resolution task which has no implementation for medical images as far as we know.

Experimental results show that the performance of CycleGAN on BraTS 2018 and on our 7T dataset is better and more stable than that of the ESRGAN model. In further works, we intend to validate our results with subjective tests made performed on a population of radiologists trained for 7T analysis.

Figure 1 :

 1 Figure 1: Architecture of building blocks in network. (i) Residual blocks in SRGAN[3] (ii) Residual blocks without BN in[START_REF] Lim | Enhanced deep residual networks for single image super-resolution[END_REF] (iii) Dense blocks in SRDenseNet[START_REF] Tong | Image super-resolution using dense skip connections[END_REF] (iv) Residual Dense Blocks from (ii) and (iii) in[START_REF] Zhang | Residual dense network for image super-resolution[END_REF], (v) Residual-in-residual dense blocks in ESRGAN[START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF]

Figure 2 :

 2 Figure 2: Network architecture of ESRGAN generator [21].

Figure 3 :

 3 Figure 3: Network architecture of CycleGAN generator. The generator contains several RDBs with local residual learning structure to improve information flow. At the higher architecture, the global residual learning structure use multi-level features to synthesize information in LR before up-sampling.

  (a) 3T T1 MRI superresolution (b) 7T T1 MRI super-resolution

Figure 4 :

 4 Figure 4: Super-resolution on MRI images. From top to bottom: full-size and zoom image. From left to right: original image, ESRGAN image, CycleGAN image

Table 1 :

 1 shows the result of ESRGAN and CycleGAN on 3T and 7T MRI images with different types of sequences. Average value of PSNR (dB) and SSIM for scale factor ×4 on 3T and 7T MRI data Table

		Bicubic	ESRGAN	CycleGAN
		PSNR SSIM PSNR	SSIM PSNR	SSIM
	3T MRI	27.26	0.52	35.85	0.88	36.76	0.92
	7T MRI	25.45	0.51	31.85	0.59	36.79	0.91


L Vaiani

E Migliorini

E A A Cavalcanti-Adam

A E E Uva

M Fiorentino

M Gattullo

V M M Manghisi

A Boccaccio
email: antonio.boccaccio@poliba.it

Coarse-Grained Elastic Network Modelling: a Fast and Stable Numerical Tool to Characterize Mesenchymal Stem Cells Subjected to AFM Nanoindentation Measurements

Keywords: Elastic Network Model, Atomic Force Microscopy, Cell Material Characterization, Meshless Methods

Coarse-grained elastic network modelling: A fast and stable numerical tool to characterize mesenchymal stem cells subjected to AFM nanoindentation measurements

Introduction

Stem cells are undifferentiated cells capable of differentiation towards specific cell phenotypes. The control of stem cells biology through mechanical factors remains a poorly understood topic and represents one of the main objectives of mechanobiology [START_REF] Lee | Stem cell mechanobiology[END_REF]. Mechanobiology is an interdisciplinary research area that integrates engineering and biology, exploiting the best computational and experimental techniques to investigate the interactions between external mechanical stimuli and internal adaptive response of living cells [START_REF] Van Der Meulen | Why mechanobiology? A survey article[END_REF] and to understand the relationships linking forces acting on a cell and biophysical responses [START_REF] Bao | Cell and molecular mechanics of biological materials[END_REF][START_REF] Lim | Mechanical models for living cells -A review[END_REF][START_REF] Humphrey | Vascular adaptation and mechanical homeostasis at tissue, cellular, and subcellular levels[END_REF]. Computational mechanobiological models were recently developed to design biomaterials for bone tissue engineering applications [START_REF] Cantore | Recent Advances in Endocrine, Metabolic and Immune Disorders: Mesenchymal Stem Cells (MSCs) and Engineered Scaffolds[END_REF][START_REF] Boccaccio | Rhombicuboctahedron unit cell based scaffolds for bone regeneration: geometry optimization with a mechanobiologydriven algorithm[END_REF][START_REF] Rodríguez-Montaño | Irregular Load Adapted Scaffold Optimization: A Computational Framework Based on Mechanobiological Criteria[END_REF].

The starting point to investigate stem cell mechanobiology is the assessment of the mechanical properties and deformation characteristics of the cell and its subcellular components. One of the most commonly adopted techniques to measure these properties is the Atomic Force Microscopy (AFM) that is widely used in experimental biology for mechanical characterization of whole living cells [START_REF] Rosenbluth | Force microscopy of nonadherent cells: A comparison of leukemia cell deformability[END_REF][START_REF] Murakoshi | Local mechanical properties of mouse outer hair cells: Atomic force microscopic study[END_REF][START_REF] Boccaccio | Nanoindentation characterisation of human colorectal cancer cells considering cell geometry, surface roughness and hyperelastic constitutive behaviour[END_REF] as well as for subcellular components [START_REF] Boccaccio | Nanoscale characterization of the biomechanical hardening of bovine zona pellucida[END_REF][START_REF] Boccaccio | A hybrid characterization framework to determine the visco-hyperelastic properties of a porcine zona pellucida[END_REF], and also for biological tissues in general [START_REF] Babu | Mechanics of brain tissues studied by atomic force microscopy: A perspective[END_REF][START_REF] Duan | How cell culture conditions affect the microstructure and nanomechanical properties of extracellular matrix formed by immortalized human mesenchymal stem cells: An experimental and modelling study[END_REF][START_REF] Herrera | From macroscopic mechanics to cell-effective stiffness within highly aligned macroporous collagen scaffolds[END_REF][START_REF] Piétrement | AFM mechanical mapping at the interface between a bioactive glass coating and bone[END_REF][START_REF] Rahmoun | Experimental characterization and micromechanical modeling of the elastic response of the human humerus under bending impact[END_REF].

It essentially consists in moving a spherical, conical or pyramidal nanoindenter fixed at the tip of a flexible cantilever of known stiffness into the surface of the material to investigate [START_REF] Vargas-Pinto | The effect of the endothelial cell cortex on atomic force microscopy measurements[END_REF][START_REF] Liu | Deformation behaviour of soft particles: A review[END_REF]. The tip displacement is detected and, based on the cantilever stiffness, is converted into force, in order to obtain a force-indentation curve that is strictly related to the stiffness of the material under examination. The effects of the AFM tip geometry and of the stress state acting on the nanoindented material were investigated in previous studies [START_REF] Boccaccio | Effect of AFM probe geometry on visco-hyperelastic characterization of soft materials[END_REF][START_REF] Boccaccio | Effect of the residual stress on soft sample nanoindentation[END_REF]. A common practice for assessing elastic properties of a biomaterial starting from force-indentation AFM curves exploits the hertzian contact theory [START_REF] Liu | Deformation behaviour of soft particles: A review[END_REF][START_REF] Lin | Spherical indentation of soft matter beyond the Hertzian regime: Numerical and experimental validation of hyperelastic models[END_REF]. The main hypotheses of this theory, in the particular case of a spherical indenter, are the following: (i) the strains in the nanoindented material are infinitesimal, (ii) the contact area is small, (iii) the nanoindented material has a linear elastic and uniform behavior, (iv) the nanoindented material presents a perfectly flat surface with infinite dimensions, (v) the indenter is perfectly spherical. In a real AFM tip to living cell contact not all those conditions are satisfied and is impossible to distinguish, by using this theoretical framework, the single contributions to the overall stiffness of the several subcellular components [START_REF] Boccaccio | Nanoindentation characterisation of human colorectal cancer cells considering cell geometry, surface roughness and hyperelastic constitutive behaviour[END_REF][START_REF] Migliorini | Nanoindentation of Mesenchymal Stem Cells using Atomic Force Microscopy: Effect of Adhesive Cell-Substrate Structures[END_REF], acting as a series-system composed by different stiffness values, but only retrieve an overall cell stiffness as a weighted average of the subcellular components stiffnesses.

Nowadays computer simulations represent a powerful tool for studying physical phenomena presenting real conditions that are not reproducible in a laboratory [START_REF] Steinhauser | A review of computational methods in materials science: Examples from shock-wave and polymer physics[END_REF]. A common numerical technique suitable for modelling deformable objects subjected to loads and constraints is the Finite Element Method (FEM) [START_REF] Huiskes | A survey of finite element analysis in orthopedic biomechanics: the first decade[END_REF][START_REF] Zienkiewicz | The finite element method in engineering science[END_REF]. The main issues encountered in the finite element analysis of soft biological materials are the following: (i) the numerical convergence is problematic and often missing even for small imposed displacements due to very low Young's moduli; (ii) the large local distortions induced in the mesh and the consequent need for numerical stabilization functions; (iii) the high computational resources required for modelling complex structures. The main goal of this study is to overcome the problems related to the large deformations and distortions imposed on mesh-based models in simulations of soft biological materials. A meshless approach is therefore proposed, characterized by much more compliant topologies, that have no continuity constraints on displacement fields and related deformations [START_REF] Rajagopal | Computational modeling of single-cell mechanics and cytoskeletal mechanobiology[END_REF]. Some examples of meshless modelling of living cells are presented in the review paper by [START_REF] Ye | Particle-based simulations of red blood cells-A review[END_REF] [START_REF] Ye | Particle-based simulations of red blood cells-A review[END_REF] where several particle-based simulations of red blood cells interacting with a viscous fluid flow are compared. Very interesting is also the work of [START_REF] Vassaux | Stem cell mechanical behaviour modelling: substrate's curvature influence during adhesion[END_REF] [START_REF] Vassaux | Stem cell mechanical behaviour modelling: substrate's curvature influence during adhesion[END_REF], who studied the effect of the curvature of substrate during the adhesion of stem cells through a particle-based approach.

In this study, AFM nanoindentation measurements were first conducted on human mesenchymal stem cells. Then, a Coarse-Grained Elastic Network Modelling (CG-ENM) [START_REF] Ingólfsson | The power of coarse graining in biomolecular simulations[END_REF][START_REF] Togashi | Coarse-grained protein dynamics studies using elastic network models[END_REF][START_REF] Kot | Elastic moduli of simple mass spring models[END_REF] was used to represent the stem cell and its subcellular components and to simulate the nanoindentation process. The pre-tensioning stress state acting on the cell due to the action of stress fibers that anchor the cell to the substrate, was also simulated and the Young's moduli of different subcellular components, such as cytoskeleton and cortex, were extracted. The results obtained with the elastic network model were then compared to those computed with the Hertz contact theory and the Finite Element Method (FEM), showing a good agreement. We found that the proposed approach generally leads to a stable numerical modelling even for large deformations and to significantly smaller computational costs than those typically required in mesh-based numerical tools.

Materials and methods

Study overview

Aiming at extracting a set of global Young's moduli, a physical-mathematical scheme has been proposed (Fig. 1), implementing a reversed approach with respect to the classical characterization of new materials performed through tensile tests. In detail, the workflow includes the following steps and procedures (Fig. 1):

(i) nanoindentation experiments;

(ii) generation of a spring network model simulating the nanoindentation process of a stem cell and identification of the optimal local material properties through numerical fitting of real forceindentation curves;

(iii) transferring of the topology and of the optimized local material properties of the spring network model to a "monodimensional twin" model;

(iv) extraction of the global elastic parameters through a second computer experiment simulating a common tensile test.

The last step allowed the comparison of the obtained global Young's moduli to the tabular data available in the literature. Each block illustrated in Fig. 1 was extensively described in the following sections.

AFM experiments

For AFM nanoindentation measurements, mesenchymal stem cells (MSCs) from human primary material were used, donated by the Professor Ho group from the Department of Medicine V, Heidelberg University and isolated as described in a previous study [START_REF] Wagner | Replicative senescence of mesenchymal stem cells: A continuous and organized process[END_REF].

After seeding of the mononuclear cell fraction evolving colonies were separated and MSC further expanded in plastic culture flasks. Cells were cultured in a sub-confluent monolayer in growth medium consisting of DMEM-LG supplemented with MCDB201 (40% (v/v), Sigma), L-glutamine (2 mM, Sigma), penicillin/streptomycin (100 U/mL, Lonza), insulin transferrin selenium (1% (v/v), Sigma), linoleic acid albumin from bovine serum albumin (1% (v/v), Sigma), dexamethasone (10 nm, Sigma), L-ascorbic acid-2-phosphate (0.1 mM, Sigma), PDGF subunit B (PDGF-BB), epidermal growth factor (EGF) (both 10 ng/mL; PreproTech, Rocky Hill, NJ, USA) and FCS (2% (v/v), HyClone, GE Healthcare). Medium was changed twice per week and early passages 3-5 were used.

MSCs were detached with Accutase (Sigma Aldrich) from the flask and plated on glass cover slips (24×24 mm 2 ; Menzel Gläser, Braunschweig, Germany), previously coated with 1 mg/mL of The stiffness of MSCs was tested one week after osteogenic induction via osteogenic media.

Indentation curves were obtained along the cell body crossing the cellular soma. The AFM was set to have a maximal indentation value no larger than 200 nm to obtain forces applied by the cantilever to the sample always less than 5 nN. For the calculation of the cellular mechanical properties, we averaged five curves taken at the highest point of the cell (corresponding to the soma) and we repeated the measurements on ten different mesenchymal stem cells. In detail, before computing the average value, the five curves were aligned in correspondence of the point where the AFM tip comes in contact with the cell. The function interp1 available in Matlab was utilized to determine, for each experimental curve, the value of indentation force at specific indentation depths (query points). Then, the average force-indentation curve was determined, for each of the ten cells, by simply computing the average value of force in correspondence of these indentation depths and was implemented in the optimization algorithm described below in the Section 2.7.

Geometrical modelling

Stem cells belong to the category of eukaryotic cells and are characterized by a complex internal structure. The subcellular components that were hypothesized to be included in the stem cell model are: the nucleus, the cytoskeleton, a network of interconnected filaments [START_REF] Blanchoin | Actin dynamics, architecture, and mechanics in cell motility[END_REF], and the cell cortex. Stress fibers and focal adhesions, that guarantee the adhesion of the cell to the substrate where it lies, were also modelled. Focal adhesions are transmembrane receptors that establish the adhesion between the ventral side of the cell and the extracellular matrix [START_REF] Oria | Force loading explains spatial sensing of ligands by cells[END_REF][START_REF] Liu | Nanoparticle tension probes patterned at the nanoscale: Impact of integrin clustering on force transmission[END_REF][START_REF] De Beer | Force-induced destabilization of focal adhesions at defined integrin spacings on nanostructured surfaces[END_REF]. Stress fibers are contractile bundles of actomyosin subjected to tensile stress, mechanically interacting with focal adhesions and additional cellular structures [START_REF] Livne | The inner workings of stress fibers -From contractile machinery to focal adhesions and back[END_REF] such as the neighboring stress fibers and the cell cortexand exerting on them a traction force of 10 nN [START_REF] Deguchi | Tensile properties of single stress fibers isolated from cultured vascular smooth muscle cells[END_REF][START_REF] Théry | Cell distribution of stress fibres in response to the geometry of the adhesive environment[END_REF]. More in detail, three main types of stress fibers exist, distributed in different regions and characterized by different lengths, from 10 to 100 µm: dorsal stress fibers, ventral stress fibers and transverse arcs [START_REF] Burridge | Focal adhesions, stress fibers and mechanical tension[END_REF].

The CAD model of the cell was developed in the Rhinoceros CAD software environment (Version 6, McNeel) integrated with the Grasshopper plugin. The main geometry is essentially a revolved surface (Fig. 2), with an overall size of 50 µm for diameter and 13 µm for height (Fig. 2c), which is consistent with the typical size of mesenchymal stem cells [START_REF] Su | Nanomorphological and mechanical reconstruction of mesenchymal stem cells during early apoptosis detected by atomic force microscopy[END_REF]. In order to reduce the computational effort and increase the resolution at the same time, only a quarter of the whole cell was modelled. Following previous studies [START_REF] Vassaux | Stem cell mechanical behaviour modelling: substrate's curvature influence during adhesion[END_REF][START_REF] Rogge | Automated segmentation and quantification of actin stress fibres undergoing experimentally induced changes[END_REF][START_REF] Nagayama | Estimation of single stress fiber stiffness in cultured aortic smooth muscle cells under relaxed and contracted states: Its relation to dynamic rearrangement of stress fibers[END_REF], nine stress fibers were modelled in the quarter model: three arcs, three dorsal (i.e. two entire and one comprised of two halves) and three ventral stress fibers. This quarter model includes the following components (Fig. 2b):

(i) the nucleus, represented by a sphere quarter of radius 5 µm, with the center placed at 6 µm from the bottom surface;

(ii) 8 tensioning regions for dorsal stress fibers, 4 located in the medium-upper part of the cell cortex (representing the cell regions where the dorsal stress fibers apply traction forces in downward direction) and 4 at the bottom of the cell (representing the regions where the dorsal stress fibers are connected to the focal adhesions) (Fig. 2b);

(iii) 4 tensioning regions for transverse arcs, located in the lower part of the cell cortex and connected to the dorsal stress fibers (Fig. 2b);

(iv) 6 tensioning regions at the bottom of the cell where the ventral stress fibers are connected to the focal adhesions (Fig. 2b);

(v) the cortex region that was separated from the cytoskeleton, thus obtaining a 200 nm thick membrane enclosing the whole cell (Fig. 2c). The value of the cell cortex thickness used in the model is consistent with that hypothesized in previous studies [START_REF] Boccaccio | Nanoindentation characterisation of human colorectal cancer cells considering cell geometry, surface roughness and hyperelastic constitutive behaviour[END_REF][START_REF] Vargas-Pinto | The effect of the endothelial cell cortex on atomic force microscopy measurements[END_REF];

(vi) the cytoskeleton, resulting as the volume obtained by subtracting, from the initial volume of the whole cell, the previous volumes.

tensioning regions for dorsal stress fibers 4 tensioning regions for transverse arcs 6 tensioning regions for ventral stress fibers

The strategy of using tensioning regions allowed to distribute the traction force (which is a relatively intense force, having the same order of magnitude as the force exerted by the AFM tip on the whole cell) exerted by stress fibers on large-dimensioned regions, thus avoiding numerical issues.

In fact, the application of the whole traction force on a single point would produce a deformation so large that any computational tool (including finite element method and coarse-grained elastic network model) would inevitably undergo non-convergence problems. Tensioning regions not only allow solving numerical issues, but also reproduce in a reasonably correct way the physics and the biology of the cell, where internal loads do not act on point-like regions with limited dimensions but on relatively wide regions, whose actual extensions are proportional to the force exerted to them [START_REF] Schwarz | Calculation of forces at focal adhesions from elastic substrate data: The effect of localized force and the need for regularization[END_REF][START_REF] Balaban | Force and focal adhesion assembly: a close relationship studied using elastic micropatterned substrates[END_REF][START_REF] Burridge | The tension mounts: Stress fibers as force-generating mechanotransducers[END_REF].

Spring network generation

The mathematical/physical approach used to model the stem cell, treated as a deformable object, is a coarse-grained elastic network modelling [START_REF] Ingólfsson | The power of coarse graining in biomolecular simulations[END_REF]. It consists of a uniform network of interconnected springs characterized by a single stiffness magnitude for the elements inside every single functional region, implementing an internal force field acting between pairs of coarse-grained lumped masses representing groups of atoms or molecules. The mutual connections between the different spring elements are spherical joints, thus enabling the maximum structural compliance, peculiar of a soft biological material.

To implement this framework, the CAD model (Fig. 3a) was exported from Rhinoceros and given in input to the open source meshing software Gmsh (version 4.5.2). A tetrahedral volumetric mesh was then generated (Fig. 3b), appropriately selecting a refined Delaunay grid algorithm as a discretization rule, aiming at maximizing uniformity of spring lengths and isotropy in behavior. The average element size was set equal to 0.350 µm which is a sufficiently small and adequate dimension to satisfactorily solve the problem of variable contact between the nanoindenter and the upper region of the cell cortex. The next operation was the transformation of the tetrahedral mesh into a spring network model (Fig. 3c). To do this, we converted the edges of the tetrahedrons generated in the meshing phase into springs and the vertices of the tetrahedrons into nodes where springs are interconnected. The information associated to nodes position and spring connectivity was processed by an ad hoc algorithm coded in the Grasshopper plugin, for the generation of the spring network into the CAD environment. The result of the conversion process is a pattern of springs arranged in a tetrahedral structure according to a refined Delaunay topology (Fig. 3c). The total number of springs was 887605 and the average rest length of springs was 0.345 µm.

To generate the quarter model, a portion of the whole cell larger than the quarter portion was initially considered (Figs. 3a-c). This portion was meshed and converted into the spring network model. A clipping procedure was finally executed (Fig. 3d) to isolate the portion of the spring network exactly corresponding to the quarter of the cell. With this strategy, we solved two main problems that arise when the quarter model is built directly using cutting planes exactly corresponding to the symmetry planes of the cell. The first problem is related to the isotropy of the model: the refined Delaunay algorithm, selected to build the tetrahedral mesh, populates the outer planar surfaces of the geometrical model with a quite regular triangular tessellation and so the density of nodes on these surfaces will be higher than that of any other plane crossing the model. This implies that if the outer planes of the geometrical model correspond to the symmetry planes, on those planes a larger density of nodes and hence a larger density of springs will result, and also this local distribution of springs will be perfectly planar, which contrasts the hypothesis of isotropy of the spring network. The second problem is related to the structural response of the nanoindented material. If the outer planes of the geometrical model correspond to the symmetry planes of the cell, at the intersection between the two planes, the refined Delaunay algorithm will build tetrahedrons with one of the edges aligned to the central axis around which the revolution of the cell model takes place. These edges will be converted into springs and due to their alignment with the central axis, which is also the direction along which the nanoindenter acts, the structural response of the cell model to the nanoindentation process will mostly depend on those aligned springs. However, this contrasts the hypothesis of isotropic mechanical behavior of the model. The above-mentioned problems are overcome and solved with the described clipping procedures that allow, definitely, leaving a quarter of the cell with an isotropic and a more "natural and realistic" springs distribution, especially in the regions close to the symmetry planes.

Definition of spring groups

The next phase was the definition, into the CAD environment, of different spring groups functional to the subsequent simulation (Fig. 4). The spring groups that were defined are the following. (i) Nucleusthe central spherical region characterized by higher stiffness compared with that of the cytoskeleton. (ii) Tensioning groupsthe regions of the cytoskeleton pulled downwards or upwards (in case of dorsal stress fibers) or pulled together (in case of transverse arcs). (iii) Cortex the external surface surrounding the whole cell, characterized by a constant thickness. (iv) Basethe lowest layer of the cell, representing the flat surface adherent to a substrate. (v) Contactthe portion of the cell cortex that enters in contact with the spherical indenter. (vi) Freethe remaining unconstrained region mainly corresponding to the cytoskeleton, affected by the boundary conditions imposed by adjacent groups, always interconnected to it.

Simulation

All the geometrical entities composing the spring network model of the cell were then imported in Abaqus (version 6.14, Dassault Systèmes) where a spherical nanoindenter with a 5 µm radiuswhich is the same radius of the nanoindenter used in the experimentswas created and modelled as a rigid part aligned to the revolution axis of the cell model.

The springs were modelled as truss elements, i.e. two-node elastic rods with a constant cross section A, set equal to 1 µm 2 , and a rest length L0 equal to the average dimension of the edges of the

tensioning groups

only along the axis connecting the nodes. As the truss elements are connected one each other with spherical joints, that constraint only the relative positions, leaving free the relative rotations, no bending stiffness has been defined for this kind of elements. The extensional spring stiffness 𝐸 𝑙𝑜𝑐 𝐴/𝐿 0 was defined by setting a local Young's modulus, one for cytoskeleton springs Eloc,cyt and one for cortex springs Eloc,cor, assuming for the nucleus region a local modulus Eloc,nuc one order of magnitude higher than the average stiffness of the cytoskeleton. Eloc,cor and Eloc,cyt were determined via the optimization algorithm outlined in the next section, that iteratively compares the numerical forceindentation curve computed in the simulation with the same curve obtained in the nanoindentation measurements previously described. The choice of setting the value of the material properties for the nucleus and not to optimize it derives from the fact that, being the nucleus rather distant from the nanoindenter tip, the system is not sufficiently sensitive to "capture" the structural response of the nucleus and therefore to analyze its mechanical behavior. The same hypothesis was followed in a previous study [START_REF] Boccaccio | Nanoindentation characterisation of human colorectal cancer cells considering cell geometry, surface roughness and hyperelastic constitutive behaviour[END_REF].

The contact between the nanoindenter and the upper region of the cell was assumed as frictionless and was solved by implementing the "hard contact" algorithm available in Abaqus.

The simulation included two different steps: tensioning of the stress fibers and indentation, consisting in the rigid translation of the spherical indenter into the cell surface. Both steps were treated on a physical-mathematical level as incremental static problems, thus excluding all dynamic contributions (inertia, time-dependent effects, dissipations, etc.). The condition of "nonlinear geometry" was enabled in every simulation, in order to correctly map the large local deformation field in the contact region.

The boundary and loading conditions acting on the model were the following (Fig. 5). (i) No displacement or rotation allowed for the base layer of the cell (encastre) for reproducing adhesion.

(ii) 10 nN [START_REF] Deguchi | Tensile properties of single stress fibers isolated from cultured vascular smooth muscle cells[END_REF][START_REF] Théry | Cell distribution of stress fibres in response to the geometry of the adhesive environment[END_REF][START_REF] Schwarz | Calculation of forces at focal adhesions from elastic substrate data: The effect of localized force and the need for regularization[END_REF] forces imposed on every single tensioning group for stress fibers (for the tensioning groups crossing the symmetry planes a 5 nN force was applied in order to achieve a 10 nN force when mirroring the quarter model). The overall force in every group was distributed as 25% in an outer shell and 75% in an inner half-sphere, in order to realize a sort of force gradient when moving from the highly-tensioned central region to the unloaded free elements neighboring the tensioning groups [START_REF] Schwarz | Calculation of forces at focal adhesions from elastic substrate data: The effect of localized force and the need for regularization[END_REF]. In order to implement such a condition, the force per node in the inner half-sphere region F1,inn was set equal to 2F1,out, the latter being the force per node for the nodes in the outer shell region (Fig. 6). With this strategy, a non-natural and non-physical sudden step-shaped magnitude change of the force acting on the tensioning groups was avoided, thus realizing a gradient for the force per node [START_REF] Schwarz | Calculation of forces at focal adhesions from elastic substrate data: The effect of localized force and the need for regularization[END_REF][START_REF] Balaban | Force and focal adhesion assembly: a close relationship studied using elastic micropatterned substrates[END_REF][START_REF] Burridge | The tension mounts: Stress fibers as force-generating mechanotransducers[END_REF]. (iii) A 200 nm displacement imposed to the rigid spherical indenter in the second step. (iv)

Symmetry constraints imposed on the symmetry planes delimiting the model.

The measurement units used in the simulation were µm for linear dimensions and nN for forces; this combination has proved to be optimal to ensure the stability of the solver, thus avoiding under/overflow issues.

Optimization

An optimization algorithm was written in Matlab environment for iteratively comparing the force-indentation curves obtained as experimental AFM reports to those predicted in the numerical simulations (Fig. 7). In detail, in the optimization algorithm, the fmincon tool available in the optimization toolbox of Matlab, was implemented, that is devoted to find the minimum of constrained nonlinear multivariable functions. Starting from guess values for Eloc,cor and Eloc,cyt, defining the local spring stiffnesses for cortex and cytoskeleton respectively and perturbing them, different resulting curves were obtained as output from the numerical solver. Following a previous study [START_REF] Boccaccio | Nanoscale characterization of the biomechanical hardening of bovine zona pellucida[END_REF], for each cell sample, optimization runs were started from five different sets of material properties previously randomly generated. Thanks to this strategy and to the large range of variability chosen for material parameters, the whole search space was covered, thus increasing the probability of finding the global optimum. The algorithm iterated this cycle until the least squared difference between experimental and numerical force-indentation curves was less than a pre-defined numerical threshold. The optimal values of Eloc,cyt,optimfor cytoskeletonand of Eloc,cor,optimfor cortexwere then given in input to

Extraction of mechanical properties from the monodimensional twin model

The same procedure previously described to build the model of the cell was replicated for a "monodimensional twin" model, i. Successively, a clipping procedure was utilized to isolate the portion of the spring network exactly delimited by the symmetry planes (Fig. 8). With this strategy, we guaranteed the isotropic structural response of the monodimensional twin model.

For each of the ten stem cell samples, two optimal values were determined through the optimization algorithm above described (sect. 2.7): Eloc,cyt,optim, for the cytoskeleton and Eloc,cor,optim

H 1DT =70 μm D 1DT =14 μm E cor 𝑜𝑟 E cyt = 4 • 𝐹 1 • 𝐻 1𝐷𝑇 𝜋 • 𝐷 1𝐷𝑇 2 • 𝛿 1 𝛿 5 =
for the cortex. These quantities were then used for defining the stiffness of the springs constituting the monodimensional twin model and hence a virtual tensile test to extract a typical forcedisplacement curve was performed. In the simulations, constraints preventing all the displacements were imposed to the nodes lying on a flat base of the cylinder, while a displacement of δ5 = 5 µm was imposed to the nodes lying on the other flat base via a rigid plate that was utilized to easily retrieve the value of the reaction force. Therefore, for each cell sample, two force-displacement curves were traced, one was obtained by assigning to the springs the Young's modulus Eloc,cyt,optim (for the cytoskeleton), the other by assigning the Young's modulus Eloc,cor,optim (for the cortex). As in the common experimental practice, the slope of the curve was extracted for the calculation of the resulting global Young's moduli Ecor and Ecyt, for cortex and for cytoskeleton, respectively. The value of force F1 registered for the displacement δ1 = 1 µm was utilized to determine the global Young's modulus as (Fig. 8).

E cor (𝑜𝑟

E cyt ) = 4•𝐹 1 •𝐻 1𝐷𝑇 𝜋•𝐷 1𝐷𝑇 2 •𝛿 1 (1)

Results and Discussion

The average experimental force-indentation curves were compared to those predicted numerically (as described in Sect. 2.7) and the optimization algorithm changed many times the elastic parameters of the spring network model of the stem cell, until the difference between the two curves was minimized (Table 1). The rather high values of the correlation coefficient R 2 , used as a metric for the detail level of the curve fitting, could be considered as a proof of the reasonably good adequacy of the model to reproduce the physics of the problem (Fig. 9, Table 2). Error bars were also superimposed to the average experimental force-indentation curves to show the dispersion of data at different indentation depths (Fig. 9).

Through the procedure described in Sect. 2.8, the global Young's moduli for cortex and cytoskeleton were then extracted from the monodimensional twin models and compared to those calculated in a previous study [START_REF] Migliorini | Nanoindentation of Mesenchymal Stem Cells using Atomic Force Microscopy: Effect of Adhesive Cell-Substrate Structures[END_REF] implementing both a FEM model based on the Arruda-Boyce hyperelastic constitutive law [START_REF] Lin | Spherical indentation of soft matter beyond the Hertzian regime: Numerical and experimental validation of hyperelastic models[END_REF][START_REF] Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF] and an analytical model based on the Hertz contact theory [START_REF] Harding | The elastic stresses produced by the indentation of the plane surface of a semi-infinite elastic solid by a rigid punch[END_REF] (Table 3). In detail, the values of Ecyt and Ecor predicted through the Finite Element Method and listed in Table 3, were computed starting from the optimal hyperelastic parameters proper of the Arruda-Boyce law, i.e. the shear moduli µ8chain,cyt (for cytoskeleton), and µ8chain,cor (for cortex) according to the following relationships: Ecor = 2(1+)µ8chain,cor, Ecyt = 2(1+)µ8chain,cyt where , the Poisson's ratio, was set equal to 0.4999 to account for material incompressibility [START_REF] Boccaccio | Nanoscale characterization of the biomechanical hardening of bovine zona pellucida[END_REF][START_REF] Vargas-Pinto | The effect of the endothelial cell cortex on atomic force microscopy measurements[END_REF].

The proposed coarse-grained elastic network model has some limitations. First, the Poisson's ratio is a non-controllable parameter [START_REF] Kot | Elastic moduli of simple mass spring models[END_REF], which directly depends on the discretization rule of the deformable volume, consisting in a nearly-regular and hence isotropic network of interconnected tetrahedrons. Second, the organelles residing into the cytoplasm were considered as a part of it and their material properties were "confused" with those of cytoskeleton. Third, even if the proposed model is capable of distinguishing, unlike the hertzian contact theory, the material properties of different subcellular components, in the current study a clear and direct experimental demonstration of the correctness of the model predictions is missing. Ideally, to demonstrate this, nanoindentation measurements should first be carried out on cells, the proposed coarse-grained elastic network model should then be implemented to determine the material properties of different subcellular components.

At this point, the different subcellular components should be physically isolated, nanoindented and analyzed, one by one, with the proposed model. If the model works correctly, the material properties predicted before and after the isolation should be the same. An intriguing procedure recently proposed [START_REF] Fischer | Effect of Nuclear Stiffness on Cell Mechanics and Migration of Human Breast Cancer Cells[END_REF] that can be adopted to separate the single contributions may consist in treating pharmacologically and selectively the single subcellular component. One should first nanoindent the cell and retrieve the force-indentation curve. Based on this curve, the proposed model should be implemented to determine the mechanical properties of the different subcellular components. At this point, a pharmacological treatment can be made to alter the mechanical properties of a specific subcellular component. Obviously, this leads to the alteration of the "global" force-indentation curve. We expect that by implementing the proposed framework to the altered force-indentation curve, the system will predict altered mechanical properties for the pharmacologically targeted component and the same mechanical properties -as those predicted in the cell before the pharmacological treatment -for the components not targeted. However, due to compensatory effect inside the cell, we believe that the selectivity of the pharmacological treatment could be an issue to investigate. For instance, the chromatin condensation in the nucleus with Trichostatin A (TSA) has been targeted in a previous work but that induced an unpredictable change in the stiffness of the nucleus that is cell-dependent and a slight stiffening of the cytoplasm as compensatory effect [START_REF] Fischer | Effect of Nuclear Stiffness on Cell Mechanics and Migration of Human Breast Cancer Cells[END_REF]. In the future, it will be interesting to address these aspects by engineering a synthetic eukaryote cell with a bottom-up approach and therefore isolating the influence of each component on the mechanical properties of the engineered cell, as a result of a multi-variable experimental plan. Fourth, hypothesizing to schematize the cell as a revolved surface, just one-quarter model was considered and symmetry constraints were applied on the symmetry planes. Indeed, symmetry constraints can give edge effects especially in the case where the extension of the constrained region is comparable to the extension of the unconstrained region.

However, in the proposed model, the extension of the constrained region is much smaller than that of the unconstrained region. Just to give an idea of this, we can state that the number of nodes where symmetry constraints act is the 5.97% of the total number of nodes. Preliminary analyses revealed that the generic force-indentation curve computed with the one-quarter model is practically overlapped to that predicted with the entire model. For instance, implementing Ecyt and Ecor computed for the cell sample 9 we found, for a 200 nm indentation depth, that the force predicted with the entire model differ from that predicted with the quarter model by less than 0.6%. This consideration allows us to conclude that the edge effects related to the use of symmetry constraints are negligible and hence the predictions of the one-quarter model are overlapping to those obtained with the entire model.

However, in spite of these limitations, the predicted displacement fields appear consistent with the nanoindenter geometry and with the physics of the problem (Fig. 10 (a)) as well as the predictions of the proposed model are consistent with the experimental results reported in the literature. The predicted material properties, in fact, are in good agreement with those measured via atomic force microscopy and Hertz contact theory in other studies [START_REF] Chen | AFM studies of cellular mechanics during osteogenic differentiation of human amniotic fluid-derived stem cells[END_REF][START_REF] Kiss | Elasticity of human embryonic stem cells as determined by atomic force microscopy[END_REF][START_REF] Liang | Probing stem cell differentiation using atomic force microscopy[END_REF]. Interestingly, the material properties predicted through the coarse-grained elastic network model are consistent with those computed via the Finite Element Method and via the Hertz contact theory reported in a previous study [START_REF] Migliorini | Nanoindentation of Mesenchymal Stem Cells using Atomic Force Microscopy: Effect of Adhesive Cell-Substrate Structures[END_REF]. It is also interesting to observe, and is a further proof of the correctness of the model predictions, that the values of the Young's modulus EHertz computed with the Hertz contact theory fallsfor each of the ten cell samples investigated and for both, coarse-grained elastic network model and FEMwithin the interval [Ecyt; Ecor], and are closer to Ecyt. This result is consistent with the physics of the problem.

Given the prevalence of the cytoskeleton volume with respect to the cortex volume, it is reasonably correct that the average material properties predicted by the coarse-grained elastic network model are closer to those of the most prevalent material. In order to make a coherent statistical comparison between the sets of Young's moduli obtained with the coarse-grained elastic network model and with the finite element method, to the values computed via the Hertz contact theory, an equivalent Young's modulus must be defined for the former cases. To this purpose, the following semi-quantitative procedure is therefore proposed. The action of the AFM tip on the nanoindented cell can be schematized, according to a previous study [START_REF] Ogar | Contact characteristics of a sphere with a layered elastic half-space[END_REF], with a force acting on two springs disposed in series (Figure 10(b)). Other semi-empirical approaches have been recently proposed to describe the mechanics of the double layer on which the nanoindenter acts [START_REF] Lüchtefeld | Elasticity spectra as a tool to investigate actin cortex mechanics[END_REF][START_REF] Doss | Quantitative mechanical analysis of indentations on layered, soft elastic materials[END_REF][START_REF] Rheinlaender | Cortical cell stiffness is independent of substrate mechanics[END_REF]. The first spring, having a length Lcor = 200 nm (i.e. the thickness of the cell cortex), a cross section Acor and an elastic constant

Kcor = E cor •𝐴 𝑐𝑜𝑟 𝐿 𝑐𝑜𝑟
, is representative of the elastic response of the cell cortex; the second spring, having a length Lcyt = 600 nm (i.e. the height of the cytoskeletal region where non-zero stress values act during the 200 nm nanoindentation), a cross section Acyt and an elastic constant Kcyt =

E cyt •𝐴 𝑐𝑦𝑡 𝐿 𝑐𝑦𝑡
, is representative of the elastic response of the cytoskeleton. Applying the classical relationships for springs disposed in series and following the hypothesis that Acor = Acyt it can be easily demonstrated that the equivalent Young's modulus Eeq of the spring that describes the elastic response of the two springs having Kcor and Kcyt as stiffnesses, can be expressed as:

E eq = (𝐿 𝑐𝑜𝑟 + 𝐿 𝑐𝑦𝑡 ) • E cor •E cyt E cyt •𝐿 𝑐𝑜𝑟 +E cor •𝐿 𝑐𝑦𝑡 (2)
Applying the Eq. ( 2), for each cell sample, and setting the values Ecor and Ecyt obtained with both, the coarse-grained elastic network model (CG-ENM) and the FEM (Table 3), we determined the equivalent Young's moduli and compared them with those computed via the Hertz contact theory (Fig. 10(c)). To make statistical inferences, we started to enquire whether data followed a normal distribution. We used the Shapiro-Wilk normality test, AS R94 algorithm, on all Young's moduli sets (i.e. the set of equivalent Young's moduli computed for coarse-grained elastic network model, the set of equivalent Young's moduli computed for FEM and the set of Young's moduli computed via the Hertz contact theory). All the three sets passed the normality test (EHertz: W=0.886, p-value=0.155;

Eeq_CG-ENM: W=0.914, p-value=0.312; Eeq_FEM: W=0.917, p-value=0.335). To evaluate the homoscedasticity, we applied the Levene test, that revealed that samples were homoscedastic (F(2, 27)= 0.467, p-value=0.632). Given these preliminary results, we performed ANOVA to compare the three Young's moduli distributions. It revealed that there is not a statistically significant difference between the distributions of the equivalent Young's moduli and the Young's moduli obtained with the Hertz theory (F(2, 27)=1.149, p-value=0.332).

For all samples, the R 2 values were close to 1 (Table 2), which means that the assumptions made on discretization topology (refined Delaunay algorithm for tetrahedrons and an adequate high spring density for the considered deformations) and on spring behavior (linear elastic) were reasonably correct for this analysis. R 2 values computed for the FEM are slightly higher than those obtained through the coarse-grained elastic network modelling. However, it is worthy to note that the correlation coefficient computed with the FEM method derives from an optimization process performed by managing four design variables, (µ8chain for cytoskeleton and cortex, and λL for cytoskeleton and cortex) while that computed with the coarse-grained elastic network model derives from an optimization process including only two design variables (Eloc,cor and Eloc,cyt). In other words, we can conclude that a 2-parameter model (coarse-grained elastic network model) provided a reasonably good and fast approximation of the results of a 4-parameter model (FEM, Arruda-Boyce).

An articulated model of the cell was introduced and described with a reasonable but arbitrary number of design hypotheses. Preliminary analyses were conducted to evaluate the impact of these hypotheses on the obtained results. For instance, the quarter model of the cell was hypothesized to include nine stress fibers. We evaluated how the cell mechanical properties change for a variable number of stress fibers: fifteen, nine, four, one, or for the case of stress fibers absent. Fig. 11 shows, for instance, the force-indentation curves obtained for the same material properties (i.e Eloc,cor,optim and Eloc,cyt,optim computed for cell sample 9, Table 1) and for different numbers of stress fibers. It can be seen that reducing the number of stress fibers leads, for the same nanoindentation depth, to increasing values of force that the cell opposes to the AFM nanoindenter. This result can be justified with the "qualitative" argument that the stress fibers create a tensile stress state on the cellular region in contact with the AFM tip that favors the tip penetration into the cell. As the number of stress fibers decreases, the entity of the tensile stress acting on the cell decreases too and consequently the nanoindenter 4 lists, for instance, Ecyt and Ecor computed for the cell sample 9.

However, the percent differencewith respect to the model including nine stress fibersnever exceeds 4.03% in the case of the cytoskeleton and 1.00% in the case of the cortex. These results lead us to conclude that, in view of a simpler description of the cell, a sort of reasonable result can be achieved even without stress fibers. However, the model should consider the effect of stress fibers especially in the case of cells including a large number of them. The nucleus stiffness was set as a constant value and not as a design variable to optimize. Its value, according to previous studies [START_REF] Boccaccio | Nanoindentation characterisation of human colorectal cancer cells considering cell geometry, surface roughness and hyperelastic constitutive behaviour[END_REF][START_REF] Deguchi | Assessment of the mechanical properties of the nucleus inside a spherical endothelial cell based on microtensile testing[END_REF] was set one order of magnitude higher than the average stiffness of the cytoskeleton. Nevertheless, recent literature has pointed out that this is not necessarily the case [START_REF] Fischer | Effect of Nuclear Stiffness on Cell Mechanics and Migration of Human Breast Cancer Cells[END_REF]. To evaluate the effect of the hypothesized nucleus material properties on the predicted values of Ecyt and Ecor, preliminary analyses were carried out with different values of Eloc,nuc. Figure 12 shows, for instance, the force-indentation curves obtained implementing Eloc,cyt,optim and Eloc,cor,optim predicted for the cell sample 9 and three different values of the nucleus material properties Eloc,nuc:

-Eloc,nuc with the same order of magnitude as the average stiffness of the cytoskeleton; -Eloc,nuc one order of magnitude higher than the average stiffness of the cytoskeleton (design hypothesis);

-Eloc,nuc two orders of magnitude higher the average stiffness of the cytoskeleton.

The force-indentation curves obtained with Eloc,nuc with one or two orders of magnitude higher than the average stiffness of cytoskeleton are practically overlapped. Larger differences can be seen instead, in the case of Eloc,nuc with the same order of magnitude as the average stiffness of the cytoskeleton. The global material properties predicted with the proposed coarse-grained elastic network model increase for decreasing values of Eloc,nuc. Table 5 lists, for instance, the values of Ecyt and Ecor predicted by the optimization algorithm for the three values of Eloc,nuc as well as the percent differencecomputed with respect to the model with Eloc,nuc one order of magnitude higher than the average stiffness of the cytoskeletonof the predicted values of Ecyt and Ecor. It can be seen that the highest value of this difference was predicted in the case of the cytoskeleton and does not exceed 5.18%. The global material properties were computed by implementing a monodimensional twin model with a cylindrical shape. It would be very interesting knowing how the utilized geometry affects the predicted material properties. To investigate this issue, two further (in addition to the main one) monodimensional twin models were built (Fig. 13(a)) with the same aspect ratio 5 : 1 and having volumes equal to 10% of the original volume. Interestingly, implementing Eloc,cyt,optim and Eloc,cor,optim predicted for the cell sample 9 and subjecting the models to a tensile test, the obtained stress/strain curves were very close one each other (Fig. 13(b) and (c)). The global material properties predicted with these two monodimensional twin models differwith respect to the original monodimensional twin modelby very small quantities never exceeding, in absolute value, 0.32% (Table 6). The most interesting feature of the proposed workflow is the intrinsic compliance of the spring network, due to the missing constraints for spring rotations at nodes. This is ensured by modelling the spring connections as internal spherical hinges, thus decoupling every pair of springs and not allowing bending moment transmission, but only axial forces. From a physical point of view, such a model can be regarded as an attempt to a micromechanical representation of weak and non-directional bonds in biological materials [START_REF] Vassaux | Stem cell mechanical behaviour modelling: substrate's curvature influence during adhesion[END_REF]. The associated mathematical description can generally consist in a series of simply linear systems in the form Ax=B for solving the static problem of a single incremental deformation, thus leading to a stable and fast formulation of the whole static problem of large deformation. Furthermore, it is not necessary to perform numerical interpolation on a continuum volumetric domain to compute expressions for displacement, deformation and other field variables as in finite element analysis. To demonstrate these features, as a proof of large deformation capability 1). Only markers are represented (and not lines) to better visualize the very small differences between the three models. in Fig. 14 (a) and (c) are presented undeformed and related deformed configurations for two monodimensional twin models built, one with the FEM method the other with a spring network. Both the models were hypothesized to be linear elastic. A flat base of the cylinder was constrained while a sufficiently large displacement was imposed to the other flat base via a rigid plate to understand until which extent the two models guarantee numerical convergence. In the FEM model, the value of the material properties implemented is Ecor determined for sample 1, whereas in the spring network model the value of Eloc,cor,optim (directly computed by the optimization process for the cell sample 1) was utilized. Interestingly, the spring network showed larger deformation than FEM until reaching the numerical non-convergence condition and thus the end of the analysis (Fig. 14b). The two curves presented the same tangent and were practically overlapped for small values of displacement but tended to diverge, due to the intrinsic different characteristics of the two models, for very large values of displacement.

Another remarkable result to point out is the comparison between the computation times necessary to complete the simulation. Generally, the duration of each analysis on a spring network model is about 1/3 of the time necessary to complete a FEM analysis with a comparable density of elements.

Conclusions

This paper describes a physical/mathematical method for extracting Young's moduli of a deformable solid structure. The study focused on nanoindented human mesenchymal stem cells and reproduced through a coarse-grained elastic network a real AFM experiment. An optimization algorithm iterated a loop until a satisfactory fit between real and simulated force-indentation curves was reached. The topological and behavioral characteristics of the simulated cells were then transferred to a monodimensional twin model and the Young's moduli for subcellular components were extracted through a virtual tensile test. Rigid plate used to easily retrieve the reaction force 1). Only markers are represented (and not lines) to better visualize the very small differences between the three models.

CG-ENM

Fig. 1 .

 1 Fig. 1. Scheme of the workflow implemented to extract the stem cells mechanical properties.

Fibronectin(

  Sigma) diluted in PBS. 20 nM of bone morphogenetic protein 2 (purchased from R&D Systems Inc. Minneapolis, MN, USA) has been added to the media to induce osteogenic differentiation of MSCs. The differentiation media was replaced once during the first week and just before the experiment. The cell indentation experiment was realized with an AFM (Nano Wizard, JPK Instruments, Berlin) combined with an inverted optical microscope (Zeiss Axiovert 200). Glass cover slips were mounted onto a stage (BiocellTM JPK) in differentiation medium, maintaining the temperature at 37 °C. Chromium/gold-coated cantilevers (CP-qp-SCONTBSG from sQUBE) with a colloidal probe (radius R = 5 µm) made of borosilicate glass were used with a nominal spring constant k between k = 0.006 -0.015 N/m, resonance frequency between 8 -13 kHz and length around 125 µm. Nanoindentation was performed by lowering the AFM tip onto the cell surface. After the contact with the cell, the AFM tip exerts an indentation force F = k•D, being D the cantilever deflection, which is registered for each position (Z) taken by the piezoelectric translator. The indentation (d) of the cell can hence be computed as d = Z-D.

Fig. 2 .

 2 Fig. 2. Exploiting the symmetry properties of the entire cell geometry (a), just one-quarter model was considered (b). (c) Principal dimensions in micrometers of the model used in the study.

Fig. 3 .

 3 Fig. 3. The CAD model (a) generated in Rhinoceros was given in input to the open source meshing software Gmsh (version 4.5.2) where a tetrahedral volumetric mesh was generated (b). The tetrahedral mesh was then converted into a spring network model (c). To generate the quarter model, a portion of the whole cell larger than the quarter portion was initially considered (a-c). A clipping procedure was finally executed (d) to isolate the portion of the spring network model exactly corresponding to the quarter of the cell. The aligned edges of tetrahedrons shown in (b) have been obtained at the intersection of the planes delimiting the model outer surface. The clipping procedure successively carried out, allows having not aligned edges (and hence not aligned springs) on the axis along which the nanoindenter acts and consequently, an isotropic behavior of the model.

Fig. 4 .

 4 Fig. 4. Spring groups defined in the model.

Fig. 5 .

 5 Fig. 5. Boundary and loading conditions imposed to the cell in the simulation environment. For the sake of clarity only the forces acting on the visible tensioning groups are indicated and not those acting on the internal ones.

Fig. 6 .

 6 Fig.6. The overall force in every group (a) was distributed as 25% in an outer shell and 75% in an inner half-sphere (b), in order to realize a sort of force gradient when moving from the highly-tensioned central region to the unloaded free elements neighboring the tensioning group.

Fig. 7 .

 7 Fig. 7. Scheme of the optimization algorithm implemented to determine the optimal local material properties for cytoskeleton and cortex.

  e. a cylinder of approximately the same volume of the cell, presenting an overall height H1DT = 70 μm, a diameter D1DT = 14 μm and characterized, according to the UNI EN 10002-1 standard (Metallic materials -Tensile Testing) by an aspect ratio of height : diameter = 5 : 1. This cylinder was utilized to extract the global Young's moduli of cytoskeleton and cortex (Fig.8). This second model presented the same geometrical (same spring network topology, cross sections and rest lengths) and structural (same local spring stiffnesses) characteristics of the cell. Also in this case, exploiting the symmetry of the problem, just one-quarter of the cylindrical volume was modelled and symmetry constraints were applied to the nodes lying on the symmetry planes. As in the cell model, outer planes of the geometrical model different from the symmetry planes were utilized to isolate a portion of the cylindrical volume larger than the quarter portion.

Fig. 8 .

 8 Fig. 8. (a) Simulation setup for the monodimensional twin model (dimensions in micrometers); (b) Magnified cross sectional view of the quarter model; c) Procedure for the extraction of global Young's moduli.

Fig. 9 .

 9 Fig. 9. Experimental force-indentation curves registered for the ten stem cell samples and related numerical curves predicted via the optimization algorithm. Error bars were superimposed to the average experimental force-indentation curves to show the dispersion of data at different indentation depths.

Fig. 10 .

 10 Fig. 10. (a) Normalized displacement field U Z /U Zmax predicted by the coarse-grained elastic network model (CG-ENM). (b) Scheme of the AFM tip-cell contact adopted to determine the equivalent Young's modulus. (c) Box plot of the equivalent Young's moduli compared with the Young's moduli computed via the Hertz contact theory.

Fig. 11 .

 11 Fig. 11. Force-indentation curves computed by the model for the same material properties (i.e. E loc,cor,optim and E loc,cyt,optim of sample 9) and for a variable number of stress fibers. The activated (in red) and the disabled (blue) stress fibers simulated in the model are shown in the figures above the diagram.

Fig. 12 .

 12 Fig.12. Force-indentation curves computed by the model for the same material properties (i.e. E loc,cor,optim and E loc,cyt,optim computed for sample 9) and for variable values of E loc,nuc .

Fig. 13 .

 13 Fig. 13. (a) Monodimensional twin models utilized to assess the effect of the geometry on the predicted material properties. Stress-strain curves obtained with the models with different geometry implementing E loc,cyt,optim (b) and E loc,cor,optim (c) computed for cell sample 9 (Table1). Only markers are represented (and not lines) to better visualize the very small differences between the three models.

Fig. 14 .

 14 Fig. 14. (a) 3D view of virtual tensile tests performed through coarse-grained elastic network model (CG-ENM) and FEM. (b) An example of forcedisplacement curves obtained with CG-ENM compared to that obtained with FEM. (c) Comparison between undeformed and deformed configurations of CG-ENM and FEM models.

  DISPLACEMENT [μm]Large displacement imposed via rigid plate U, U Z[µm]

Fig. 1 .

 1 Fig. 1. Scheme of the workflow implemented to extract the stem cells material properties.

Fig. 2 .

 2 Fig. 2. Exploiting the symmetry properties of the entire cell geometry (a), just one-quarter model was

Fig. 3 .

 3 Fig. 3. The CAD model (a) generated in Rhinoceros was given in input to the open source meshing

Fig. 4 .

 4 Fig. 4. Spring groups defined in the model.

Fig. 5 .

 5 Fig. 5. Boundary and loading conditions imposed to the cell in the simulation environment. For the

Fig. 6 .

 6 Fig.6. The overall force in every group (a) was distributed as 25% in an outer shell and 75% in an

Fig. 7 .

 7 Fig. 7. Scheme of the optimization algorithm implemented to determine the optimal local material

Fig. 8 .

 8 Fig. 8. (a) Simulation setup for the monodimensional twin model (dimensions in micrometers); (b)

Fig. 9 .

 9 Fig. 9. Experimental force-indentation curves registered for the ten stem cell samples and related

Fig. 10 .

 10 Fig. 10. (a) Normalized displacement field UZ/UZmax predicted by the coarse-grained elastic network

Fig. 11 .

 11 Fig. 11. Force-indentation curves computed by the model for the same material properties (i.e.

Fig. 12 .

 12 Fig. 12. Force-indentation curves computed by the model for the same material properties (i.e.

Fig. 13 .

 13 Fig. 13. (a) Monodimensional twin models utilized to assess the effect of the geometry on the

Fig. 14 .

 14 Fig. 14. (a) 3D view of virtual tensile tests performed through coarse-grained elastic network model

OPTIMIZATION (sect. 2.7) GEOMETRICAL MODEL (sect. 2.3) SPRING NETWORK MODEL (sect. 2.4 -2.5) SIMULATION (sect. 2.6) E loc,cor,optim E loc,cyt,optim E cor E cyt AFM EXPERIMENTS (sect. 2.2) GEOMETRICAL MODEL (sect. 2.8) SPRING NETWORK MODEL (sect. 2.8) SIMULATION (sect. 2.8)

Table 1 .

 1 Optimal local material properties predicted via the optimization algorithm.

		Coarse-Grained Elastic Network
	Sample N.	Eloc,cyt,optim [Pa] Eloc,cor,optim [Pa]
	Sample 1	160	508
	Sample 2	252	490
	Sample 3	208	504
	Sample 4	110	470
	Sample 5	95	415
	Sample 6	127	485
	Sample 7	230	480
	Sample 8	221	507
	Sample 9	140	490
	Sample 10	120	470

Table 2 .

 2 Correlation coefficient R

2

values related to the experimental-numerical curve fitting, computed for the proposed coarse-grained elastic network model versus the same coefficient values computed for the FEM model

[START_REF] Migliorini | Nanoindentation of Mesenchymal Stem Cells using Atomic Force Microscopy: Effect of Adhesive Cell-Substrate Structures[END_REF]

.

Table 3 .

 3 Global material properties computed with the proposed coarse-grained elastic network model compared to those predicted via the FEM method and the Hertz contact theory[START_REF] Migliorini | Nanoindentation of Mesenchymal Stem Cells using Atomic Force Microscopy: Effect of Adhesive Cell-Substrate Structures[END_REF] (applied to the experimental force-indentation data).

		Coarse-Grained	Finite Element Method,	Hertzian Contact
		Elastic Network	Arruda-Boyce	Theory
	Sample N.	Ecyt [Pa]	Ecor [Pa]	Ecyt [Pa]	Ecor [Pa]	EHertz [Pa]
	Sample 1	2911	9216	3519	9754	4110
	Sample 2	4578	8898	5290	9687	6350
	Sample 3	3771	9144	4462	9671	5200
	Sample 4	1999	8352	2575	9738	3195
	Sample 5	1722	7531	2258	9645	2800
	Sample 6	2286	8751	2901	9570	3410
	Sample 7	4152	8666	4638	9807	5550
	Sample 8	3988	9161	4551	9729	5680
	Sample 9	2527	8845	3129	9555	3455
	Sample 10	2166	8484	2640	9618	3050

Table 4 .

 4 Global material properties computed for the cell sample 9 and for different numbers of stress fibers.

		Optimal global properties Percent difference	Percent
				cytoskeleton	difference cortex
	# of stress fibers	Ecyt [Pa]	Ecor [Pa]	[%]
	15	2594	8759	2.65	-0.97
	9	2527	8845	-	-
	4	2497	8815	-1.18	-0.34
	1	2448	8770	-3.13	-0.85
	0	2425	8756	-4.03	-1.00

Table 5 .

 5 Global material properties computed for the cell sample 9, for different nucleus material

			properties Eloc,nuc
		Optimal global properties Percent difference	Percent
			cytoskeleton	difference cortex
	Eloc,nuc	Ecyt [Pa]	Ecor [Pa]	[%]
	One order of
	magnitude higher	2527	8845	-	-
	(design hypothesis)
	Two orders of magnitude higher	2523	8833	-0.16	-0.13
	Same order of magnitude	2658	8971	5.18	1.42

Table 6 .

 6 Global material properties computed for the cell sample 9, for different geometries of the

		monodimensional model
		Optimal global properties Percent difference	Percent
				cytoskeleton	difference cortex
	Volume	Ecyt [Pa]	Ecor [Pa]	[%]
	Original Volume V	2527	8845	-	-
	0.9×V	2519	8818	-0.32	-0.31
	1.1×V	2529	8853	0.08	0.09

Acknowledgements

The support given by Dr. Patrick Horn (Department of Medicine V, Heidelberg University), for donating the human mesenchymal stem cells, is gratefully acknowledged. We thank the European Union's Framework Program for Research and Innovation Horizon 2020 (2014-2020) under the Marie Sklodowska-Curie Grant Agreement No. 658334 for the funding received.

R 2 Sample N. Coarse-Grained Elastic Network FEM Sample 1 0,9962 0,9825 Sample 2 0,9836 0,9979 Sample 3 0,9862 0,9984 Sample 4 0,9776 0,9920 Sample 5 0,9733 0,9947 Sample 6 0,9746 0,9959 Sample 7 0,9929 0,9936 Sample 8 0,9893 0,9963 Sample 9 0,9952 0,9833 Sample 10 0,9949 0,9737

The resulting Young's moduli showed good agreement with the corresponding moduli resulting from a previous study, performed using the Finite Element Method, and also presented the same order of magnitude as the elastic moduli predicted by hertzian contact theory and related literature.

The proposed approach can be considered as a useful and fast tool to study the mechanical behavior of biological soft materials. As a general rule, in a coarse-grained elastic network model, the time required to complete the simulation is about 33% of the time required to run the same simulation with a FEM model presenting the same discretization density for the elements.

Remarkable numerical stability and computational speed make this method comparable to a traditional FEM analysis and even preferable in the case of large displacements and deformations of soft materials.

Declaration of Competing Interest

The authors have declared that no competing interests exist.

Ethics statement

Bone marrow from healthy donors for allogeneic transplantation was taken after written consent using guidelines approved by the Ethic Committee on the use of Human subjects at the Heidelberg University.