An Embodied Brain Model of the Human Foetus

Cortical learning via sensorimotor experiences evoked by bodily movements begins as early as the foetal period. However, the learning mechanisms by which sensorimotor experiences guide cortical learning remain unknown owing to technical and ethical difficulties. To bridge this gap, we present an embodied brain model of a human foetus as a coupled brain-body-environment system by integrating anatomical/physiological data. Using this model, we show how intrauterine sensorimotor experiences related to bodily movements induce specific statistical regularities in somatosensory feedback that facilitate cortical learning of body representations and subsequent visual-somatosensory integration. We also show how extrauterine sensorimotor experiences affect these processes. Our embodied brain model can provide a novel computational approach to the mechanistic understanding of cortical learning based on sensorimotor experiences mediated by complex interactions between the body, environment and nervous system.


Supplementary Methods
Musculoskeletal body model. We constructed the bone, muscle and skin of the foetus model mainly based on foetal MRI data and used other complementary data when necessary because of the limited resolution of the MRI data. The human foetus examined by MRI was a historical specimen with a gestational age of 206 days belonging to the Kyoto Collection 1 . The scan was conducted using a 1.5 T MRI system (Excelart Vantage, Toshiba Medical Systems, Tokyo, Japan). We used T1-weighted images and a 3D gradient-echo sequence with the following parameters: a time repetition/time echo of 30/7 ms, an imaging matrix of 256×192×176 pixels, and a spatial resolution of 0. (http://www.boneclones.com/) and acquired computed tomography scans with a spatial resolution of 0.02 mm. The MRI provided information on the global bone shape, including cartilage, although we could not obtain detailed structures of some features, such as the bones of the hands and skull owing to resolution limitations. In contrast, the data extracted from the foetal skeleton replica had high spatial resolution, but little cartilage. We constructed the foetal skeletal model by combining both data sets using MAYA (Autodesk). Finally, we confirmed that the bone sizes were within the normal range of human foetuses at 32 gestational weeks with respect to the following parameters: biparietal diameter, head circumference, humerus, radius, ulna, femur, tibia and fibula 2, 3 .
We modelled 830 muscles as piecewise line segments defined by attachment and relay points based on Lee's work 4 and built the foetal musculoskeletal body model by combining the foetal MRI data and an adult musculoskeletal model 4 . We extracted the attachment and relay points of 216 muscles from the foetal MRI data. We obtained those of the other muscles by manually transforming the adult bones to foetal bones using MAYA. We also calculated the muscle cross-sectional area of the 216 muscles extracted from the foetal MRI data to estimate the maximal voluntary contraction force 5 . We calculated the cross-sectional areas of the muscles that could not be extracted from the foetal MRI data by scaling the adult data using the average ratio calculated using the extracted foetal muscles. Because we subdivided the body model into 21 rigid-body parts for physical simulation, we excluded all muscles whose attachment and relay points belong to only one rigid-body and ultimately adapted 390 muscles in the musculoskeletal body model ( Fig. 1b and  A total of 20 joints with 36 degrees of freedom were modelled in the entire body by excluding the fingers and toes and simplifying the vertebrae. The range of motion of each joint was set such that it did not exceed the human neonate data for the limbs 6 or the adult data for the trunk, neck and scapula 7 , which were used because no neonatal data were available for these three factors. To simulate the three-dimensional rigid-body dynamics, we used the Open Dynamics Engine, which is a widely used open-source physics engine (http://www.ode.org/). Muscle dynamics and proprioceptive sensory feedback for the muscle spindle and Golgi tendon organ were modelled based on experimental data 8 . The dynamics in the original models were defined for one specific muscle; thus, we applied them to the muscles in the entire body by normalising the maximal force, length and range.
We allowed the foetal model to generate random movements within the range of motion of the joints and determined the range of lengths of each muscle to normalize the sensitivity of the muscle spindle model. Spinal neural circuit model. We employed the spinal neural circuit model 9 (Fig. 1g), which is based on experimental data and is scalable for the whole-body musculoskeletal model 10,11 . The spinal circuit model independently controls each muscle within an elementary circuit consisting of the following components: neural oscillator, sensory interneuron, and α and γ motor neurons. The spinal circuit relays sensory feedback from the muscle spindle and tactile mechanoreceptor models to cortical neurons via sensory interneurons and modulates the muscle activation and sensitivity of the spindle with basic neuromuscular loops, such as the stretch reflex and Alpha-Gamma linkage. The neural oscillators are considered the neural bases of spontaneous movements 12 . These spontaneous movements can be observed as early as muscles function 13 , and many animal studies suggest that their essential role in early nervous and motor development begins in the foetal period 14,15,16 . However, the spinal neural circuit that generates whole-body movements is still unknown. Therefore, we built a minimal spinal circuit model based on experimentally supported interneuronal connectivity and investigated whether the simulated movements could capture the normal features that have been reported in human studies.
The dynamics of the neural oscillator model are represented by the Bonhoffervan der Pol (BVP, or FitzHugh-Nagumo) equation as follows 9,17 : where and are the inputs from the cortex and spinal sensory interneurons, respectively. The corticospinal input is modelled as a constant value, = 0.6, because the mechanism by which cortical learning alters generated movements is outside the scope of this paper and because the functional development of corticospinal tracts, such as the transmission of motor signals, is not thought to begin during the early periods of cortical learning targeted in this paper 18 . The other constant parameters were set as = 0.7, = 0.675, = 1.75, = 0.013, = 0.022, and = 0.2 based on previous work 9 .
The sensory interneurons and α and γ motor neurons were modelled using the following transfer functions based on experimental data 8 : where was set to 1.5 p/s/nA 19 . For more details, see the original paper 9 .
Intrauterine and extrauterine environment models. In the intrauterine environmental condition, the centre of the uterus was spatially fixed and the centre of the abdomen of the foetal model was connected to the centre of the uterus by a ball joint. Thus, the abdomen body part could rotate only around its centre. In the intrauterine condition, the foetal model received forces from the uterine membrane, amniotic fluid and physical contact between body parts in addition to gravity and buoyancy forces. The forces from the uterine membrane and amniotic fluid were calculated for each point on the skin surface.
To simulate the extrauterine environmental condition, we placed the body model on a flat, rigid bed. The model was subject to only the forces of gravity and physical contact between body parts and the flat, rigid bed.
In both conditions, we set the acceleration of gravity to 9.8 m/s 2 . The forces due to the physical contact were calculated by the physical simulator.
In the intrauterine condition, we also simulated buoyancy forces using the Here, I (,E is the distance to object along a normal vector , ( . H is the minimum distance to which the physical contact force is distributed, and we set H = 5 mm. The 3,000 tactile mechanoreceptor models were distributed according to human two-point discrimination data 23 ( Fig. 1d and Supplementary Table 2).

Analysis of whole-body movements.
Regarding the chaotic and fractal properties of whole-body movements, we simulated foetal movements for 1,000 s and then analysed the trajectories of the four limbs based on experimental studies of human neonates and infants. To investigate the chaotic properties, we calculated the largest Lyapunov exponent 24 . When the exponent was positive, the trajectory exhibited characteristics of deterministic chaotic properties. We quantified the fractal properties by calculating the scaling exponent using detrended fluctuation analysis 25 . When the exponent was close to one, the trajectory exhibited fractal-type, long-range correlations. To investigate whether the generated whole-body movements were well-coordinated, a qualitative feature reported in human studies 26,27 , we analysed inter-limb coordination by measuring the degree of phase synchronization between muscles using the standard shuffle-corrected phase synchronization index 28 .
The whiskers in the box plots indicate the upper and lower quartiles.