In our project, the self-healing material will be applied to the surface of soft robot to achieve self-healing when the soft part of the robot is damaged during the work process. For this purpose, we designed a path-finding obstacle avoidance & interference algorithm to simulate the working state when the robot works underwater. We created a simulation map derived from the actual terrain, incorporating various factors to replicate the interference experienced by the underwater robot. By setting up experiment group with different interference incidence rate, we measured the damage rate for multiple times. The experiment demonstrates that working robots are prone to encounter obstacles and get impaired even in interference-free environment. So, the self-healing material we designed is believed to be obligatory. Subsequently, since TRn plays a vital role in self- healing, the hydrogen bond network is simulated by us to study how repetition times affect the stability of it. We found that the more times TRn repeats, the more stable the protein structure is, and then we confirmed our conclusion with the perspective of protein energy. Therefore, using cmRNA to repeatedly translate proteins can improve the self-healing effect of the final material. Considering that, L-DOPA is the key substance to achieve adhesion between materials and soft robots, we used the Michaelis equation to calculate the production of L-DOPA at different reaction times and conducted stability analysis, and obtained the optimal reaction time to improve experimental efficiency. Finally, in order to quantify the advantages of self-healing materials, we used the renewal reward theorem to calculate the optimal maintenance strategy and minimum cost of soft robots. As a result, we found that using self-healing materials can save costs.

Although the current satellite's remote sensing and multi-beam detection technology are very advanced and can detect the macroscopic topography of the seabed, it is slightly insufficient to detect some biological communities that are relatively microscopic. For example, coral may form coral reefs in warm shallow seas, while cold-water corals are more likely to form reefs and forests in deep seas. For example, under the sea surface of the Qingshui Cliff, a scenic spot in the north of Hualien County, Taiwan Province, China, there is a dense coral forest as shown in Fig. 1a.

Fig. 1 | Coral reefs topographic map real(Image source:a:https://www.vjshi.com/watch/15562010.html;b:https://pixabay.com/zh/photos/diver-underwater-ocean-swim-fishes-123286/;c:https://pixabay.com/zh/photos/underwater-coral-fish-sea-reef-5310424/)

It is possible that soft robots will replace divers in similarly complex terrains, such as exploring biological communities or collecting biological samples. However, it is obvious that in such terrains, soft robots may be more likely to collide with obstacles and suffer after being affected by some random factors (e.g. fish schools, ocean currents). So we designed an obstacle avoidance algorithm and generated random terrains to simulate the working state of the robot.

1. The speed of sound propagation in water is fixed;

2. The interference caused by various factors in the water (such as fish schools and ocean currents) only affect the position of the robot and do not damage the robot; position of the robot and do not damage the robot;

3. The robot will only be damaged if it hits a mountain.

Table 1 | Symbol Description

1. Set random seeds to ensure that the map is consistent in each simulation.

2. Generate a random number matrix $A_{n n}\left(a_{i j}\right)$ with a range of 0~1 and a normal distribution.

3. Flatten the terrain, set the number of loops (flatness), and set $a_{i j}$ to the average of the eight surrounding cells in each loop.

4. To generate a map of a specific terrain scale, we need to create a plane of appropriate height and assign values below the plane height to 0, so that we can get a terrain map of a specific scale.

Fig. 2 | Flowchart of Map Algorithms

We know that underwater soft robots may mainly process images and information through "sound vision" to extract features, then identify targets. In contrast, other sense signals attenuate too quickly in the water such as light and radio waves. Based on this fact, the principle of our designed algorithm is: the robot uses sonar to scan the fan-shaped area in the forward direction to detect the distance between the target and the sensor. Then the robot adjusts its forward direction according to feedback, either forward or backward or left or right, and keeps approaching the expected destination while avoiding obstacles.

TThe algorithm first inputs parameters such as the initial position and direction and initializes the speed direction and steering. Then search for obstacles within the sector area in the forward direction. If no obstacles are found, the search radius will be expanded. According to the positional relationship between the robot and the nearest obstacle point, the current speed direction is calculated, and then the forward direction and position at the next moment are obtained. Before calculating the next location, another short-range sector scan is performed to ensure that there is no collision with obstacles during the displacement process. The above process is repeated until an obstacle is encountered or the destination is reached.

Fig. 3 | Flowchart of Obstacle Avoidance Algorithm

1. Set three vectors: motion direction unit vector $V_{x, i}$, end direction unit vector $V_{z, i}$, and steering unit vector $V_{b, i}$, with coefficients $\beta$, $\delta$, and $\theta $, respectively. $V_{x \ 0}=V_{z \ 0}$.

2. Input the start point $O(x_0,y_0,z_0)$, the end point coordinates $E(x_{end},y_{end},z_{end})$, set the current robot position point ${O}^{\, \prime}=(x_0^{\, \prime},y_0^{\, \prime},z_0^{\, \prime})$, motion moment $t_{i}$, $V_{z \ 0}=\overrightarrow{OE}$, and map data $map(x,y,z)$.

3. Solving $V_{z, i}\&{O}^{\, \prime}(t_{i + 1})$

3.1. $t_{i}$ when the direction of motion $V_{x, i - 1}$.

3.2. Retrieve all the coordinate points $Temp(x,y,z\,)$ in the direction of $V_{x, i - 1}$ with $R=r_{0}+p_{0}\alpha $ ($\alpha $ as the speed of the robot; $p_{0}$ as the pre-set runnig time of the robot, mainly used to measure the detection radius and speed of the relative relationship) as a sector of the detection radius (according to the principle of the underwater soft robot detecting the underwater environment through sonar). Inside all the coordinates of $Temp(x,y,z\,)$, if there is no, then $R=R+1$, expand five hundred units is still no (indicating anomalies) then $V_{x, i}=-V_{x, i - 1}$, into 3.9.

3.3. Calculate the distance $\mathrm{d}$ from each point in the point set $Temp$ to ${O}^{\, \prime}$, $D=min\{d\leq R|Temp,{O}^{\, \prime}\}$, and if $D=\{\emptyset\}$. then there is no risk of crash,$V_{b, i}=0$, and enter 3.7.

3.4. If $r_1 \lt D \lt R$ ( $r_{1}=r_{0}+\frac{p_{0}}{5}\ \alpha$ is the maximum safe distance, below which there is a risk of a crash), $T(tempx,tempy\,)=min\{(x, y\,) | r_1 \lt D \lt R \,\} $.

3.4.1 Make $V_{b,i}\cdot \overrightarrow {O^{\, \prime}T}=-1$, and normalize to get the steering unit vector $V_{b, i}$.

3.4.2. At this point $V_{b, i}$ has two solutions$V_{b, i}\,(1)$, $V_{b, i}\,(-1)$, and if $acos(\frac{V_{x, i - 1} \cdot V_{b, i}\,(1)}{|V_{x, i - 1}|\cdot|V_{b, i}\,(1)|}) \lt \frac\pi2$, $V_{b, i}=V_{b, i}\,(1)$; and conversely, $V_{b, i}=-V_{b, i}\,(1)$.

3.5. if $r_{2} \lt D \lt r_{1} $ ($r_{2}=r_{0}+\frac{p_{0}}{10} \alpha$ is the minimum safe distance, below which there is a very high risk of a crash), make $V_{b, i}\cdot V_{x, i - 1}=- 1$, normalized to find $V_{b, i}$.

3.5.1. 80% chance it's $V_{b, i}=\{V_{b, i}|acos(\frac{V_{x, i - 1} \cdot V_{b, i}\,(1)}{|V_{x, i - 1}|\cdot|V_{b, i}\,(1)|}) \lt \frac{\pi}{2}\}$ ,20% chance it's $V_{b, i}=\{V_{b, i}|acos(\frac{V_{x, i - 1} \cdot V_{b, i}\,(1)}{|V_{x, i - 1}|\cdot|V_{b, i}\,(1)|}) \ge \frac{\pi}{2}\}$.

3.5.2. $V_{x, i}=- (V_{x, i - 1}+V_{b, i})$ , normalized to obtain the direction-of-motion unit vector $V_{x, i}$, into 3.10.

3.6. If $D \lt r_{2} $, $V_{x, i}=-V_{x, i - 1}$, enter 3.10.

3.7. $V_{z, i}=O^{\, \prime}E$.

3.8. $V_{x, i}=\theta V_{b, i}+\delta V_{z, i}+\beta V_{x, i - 1}$, Normalization yields a unit vector in the direction of motion $V_{x, i\,}$.

3.9. Retrieve all coordinate points in the $V_{x, i}$ direction within a sector of size $\alpha $ as the detection radius $Temp \, ''=\{(x, y, z)|map\}$, $d \, ''=\min \left\{ \overrightarrow{Temp \, ''O \, \prime}|k=1,2,3 \right\} $ is the distance from these points to ${O}^{\, \prime}$ if $d \, ''\le\alpha$, enter 3.6.

3.10. ${O}^{\, \prime}(t_{i})+\alpha V_{x, i}={O}^{\, \prime}(t_{i + 1})$.

4. Go back to 3 until $\mid\overrightarrow{O^ {\, '} E}\mid\leq min\_area$ or collision detection is triggered.

5. crash detection

5.1. Retrieve point $Temp^{\, \prime}=\{(x, y, z) |x_t-r_0\leq x \leq x_t+r_0, y_t-r_0 \leq y \leq y_t+r_0, map \}$ in a circular region of radius $r_0$.

5.2. Calculate the minimum distance from the midpoint of $Temp^{\, \prime}$ to ${O}^{\, \prime}$, $d^{\, \prime}=min\{\overrightarrow{Temp\, ' (k)O\, '}|,k=1,2,...\}$, if $d^{\, \prime}\leq r_{0}$, end the loop, otherwise continue that loop.

Key parameter settings

Table 2 | Key Parameter Settings

Simulation map

To verify the algorithm, we first used the random terrain scale map algorithm and set the number of seeds to randomly generate a map of size 1000m $\times$ 1000m.

Fig. 4 | Simulation map display

The results of multiple simulation experiments, conducted with the obstacle avoidance algorithm, are presented in Fig. 5. The presented results were randomly selected.

Fig. 5 | Robot trajectory map without interference

In 100 simulations, the collision detection mechanism embedded in the code was never triggered, suggesting that the obstacle avoidance algorithm devised by our team is capable of avoiding seafloor obstacles effectively, enabling the robot to reach the predefined endpoint within the randomly generated seafloor terrain map. While the algorithm may not be optimal, and the path may not be the shortest, the ability of the robot to avoid obstacles and reach the endpoint is sufficient for this simulation.

Real Map

To demonstrate the generalizability of our algorithm and the persuasiveness of our simulation, we sought to locate a scanned image of a genuine coral reef seafloor. The obstacle avoidance interference algorithm was run by extracting the black pixel points on the image and calculating their ratio to all pixel points, i.e., the shading ratio, as a proportion of the real terrain and the characteristics of this type of seafloor environment. This was done in order to generate a realistic simulated seafloor terrain map.

Fig. 6 | Topography of coral reefs on the seafloor of the northeastern Atlantic Ocean.(Wang et al. 2019)

Fig. 6a depicts the seabed coral reef topography in the northeastern Atlantic Ocean , with the most prominent features highlighted in black and white (Fig. 6b). The proportion of black and white pixels was calculated using self-written codes, resulting in a value of 51.6927%. This figure was then entered as a scale parameter into the random map code, which was used to generate a simulated seabed terrain map that reflected the actual terrain proportion.

Fig. 7 | Real scale map display

Two sets of start and end points were established from the bottom to the top and from the right to the left to assess the algorithm's applicability in different directions on the same map. The height was set to 115 meters, and the simulation was conducted 50 times. The results of the two trials were selected and are presented in Fig. 8. To demonstrate the versatility of the algorithm, we conducted a series of experiments, varying the preset heights while maintaining the starting and ending points. This approach allowed us to assess the algorithm's performance when applied to different height scenarios. Specifically, we set the heights to 95, 100, 110, and 120 meters, executing 50 sets of data for each. The corresponding heights were randomly selected for each set, and the obstacle avoidance algorithm was initiated without triggering collision detection. The outcomes of varying heights are illustrated in Fig. 9 (the colors are green, yellow, magenta, and red from bottom to top), which demonstrate that our algorithm exhibits robust performance.However, it is worth mentioning that there is a limit to the robustness of any algorithm.

Fig. 8 | Algorithm results for different orientations at 115 meters

Fig. 9 | Algorithm results for different heights

Subsequently, 15 pieces of seabed terrain were randomly selected from the real ocean on the world map and replicated 1:1 and then scaled down to 1000m $\times$ 1000m in size. The algorithm was then run separately without any crashes. This demonstrates that the algorithm is capable of assisting the robot in avoiding obstacles, irrespective of the terrain whether real or simulated in the absence of interference. However, it should be noted that simulated terrain may present a greater degree of complexity. The outcomes of executing our algorithm on a randomly selected terrain are illustrated in Fig. 10.

Fig. 10 | Demonstration of real map Algorithm operation(West longititude from 21.3926 to 17.3496 degree,north latitude minus from 40.1177 to 44.0903 degree,Our seafloor mapping data comes from the links below:https://download.gebco.net/)

Considering that many factors may affect the robot's state underwater, such as ocean current and fish interference, and so on. These factors will bring challenges to our robots' work. Because the real interference situation on the seabed is complicated and difficult to describe quantitatively. To simulate this mechanism, we simplify the problem by abstracting these complex and uncertain factors into random displacements. Each time the obstacle avoidance algorithm is run, a random interference event will occur with a certain probability (interference incidence). Different interference events will cause different degrees and directions of displacement to the current position of the robot, which may cause the robot to hit an obstacle, or continue to swim without any problem.

Fig. 11 | Flowchart of Disturbance Algorithm

1. Set variable of interference incidence $disturb$.

2. Set i levels of interference $DTX(i\ ) , i=1,2,\ldots $.

3. Set the probability of occurrence of the corresponding event $DT(i\ ) , i=1,2,\ldots $.

4. Generate three uniformly distributed random numbers over $[-1,1]$, normalized to obtain a random 3D unit vector $\text{e}$ that may point in any direction.

5. Generate a random number in the range $0\sim1$ that follows a uniform distribution and compare it with $disturb$ to determine whether the algorithm occurs or not.

6. Generate a random number in the range $0~1$ that is uniformly distributed, observe the interval that the number falls into, and determine which interference event occurs.

7. Return the interference displacement vector $Disturb=DTX(i\ )\cdot e$.

We set up six random interference events, each with its own probability of occurrence and corresponding interference level. We set the probability of occurrence of these six events $\mathrm{DT(i\, )~(i=1..6)}$ to 0.40, 0.25, 0.20, 0.10, 0.04, 0.01 respectively; and the interference level of the six events $\mathrm{DTX(i\ )~(i=1..6)}$ to 2, 4, 6, 8, 16, 32, 64 (in meters) respectively; The basic principle of parameter setting is that the higher the probability of occurrence of interference events, the lower the interference level.And vice versa.

In all subsequent experiments involving randomized interference algorithms, we will use the parameters set above.

TTo facilitate subsequent experiments, we have incorporated random interference into the obstacle avoidance algorithm, in conjunction with the randomized simulation maps described previously.

Fig. 12 | Flowchart of Obstacle Avoidance & Disturbance Algorithms

1. Starting from $t_{i}=2$, before each loop, make $O^{\, \prime}=O^{\, \prime}+Disturb$.

2. $Temp^{\, \prime}=\{(x, y, z) |x_t-r_0\leq x\leq x_t +r_0, y_t-r_0\leq y\leq y_t+r_0,map\}$.

3. $d^{\, \prime}=min\{\overrightarrow{Temp\, ' (k)O\, '}|,k=1,2,...\}$, if $d\, '\leq r_{0}$, end loop.

4. Make $c_{z0}=\frac{\left(c_{z0}+c_{zend}\right)}{2}$ so that the robot's current altitude is as close as possible to the goal altitude, and go to the next loop.

To simulate the situation of robot collision under a certain interference probability, we conducted multiple collision simulation experiments under the condition of 30% interference rate (we assume that collision may be more significant under this probability) in random terrain scale map. The red line represents the motion trajectory without interference, the blue line represents the collision trajectory under the interference probability of 30%, and the red ball represents the robot's explosion when it hits an obstacle (not necessarily explosive in reality, we exaggerate here). We randomly selected data from multiple groups of collisions as a display, as shown in Fig. 13.

Fig. 13 | Comparison of robot paths with and without interference

In order to simulate the impact of the soft robot on the robot itself after being affected by different factors in the real coral forest area on the seabed, we utilized our own developed obstacle avoidance & disturbance algorithms. The simulated seabed terrain was generated according to the aforementioned scale of the real terrain. Different groups were set up based on varying incidence of interference. Each group was subjected to 10 repetitions of the experiments, with the algorithms run 50 times in each experiment. The resulting statistics count the crashes for each experiment.

Although the robot may still work after a collision, it is often along with performance degradation, component damage, and even the need for salvage and repair. We believed that continuing observing their obstacle avoidance process in this case is not representative, so we only count the first collisions in each experiment. In addition, since the various interference factors on the seabed are too complex and changeable in reality, it is difficult for us to determine the specific interference situation in sea area. Therefore, we set different interference rates to cover seabed interference situations as much as possible, in order to illustrate the damage of the robot under different interferences

In the event that the random interference parameter remains unaltered, the interference incidence was set at 5%, 10%, 20%, and 30%, and designated respectively as A, B, C, and D, for a total of four experimental groups. These groups were utilized for experimental simulation. The results are presented below:

Table 3 | Damage rate statistics for each interference incidence

Fig. 14 | Trend chart of damage rates

As illustrated in Table 3, the mean value of the robot's impairment rate is remarkably rising in conjunction with an escalation in the interference incidence and an expansion in the uncertainty of the sea area. Consequently, the overall trajectory of the variance in its impairment rate is also exhibiting a gradual increase. Upon visualizing the data, it becomes evident that the impairment rate exhibits a pronounced increase within 10% to 30%, with a particularly notable acceleration observed within the 20% to 30% range. This phenomenon is illustrated in Fig. 14.

For calm seas (Group A), although the overall damage rate has a low mean, the variance is relatively large, and the actual damage rate may fluctuate around the mean value. For relatively calm sea areas (Groups B and C), the damage rate can be as high as 24%. Due to higher uncertainty, the variance is larger than Group A, and the actual damage rate may exceed this value. We believed that in one cycle, the salvage and repair costs caused by this damage frequency were unacceptable to robot companies. In conventional sea areas (Group D), the mean damage rate rises sharply (up to 42%), resulting in a corresponding increase in costs.

As for the case where the interference rate is greater than 30%, due to the increase in randomness, the mathematical expectation of the damage rate will inevitably continue to rise.

The above results show that it is very common for underwater robots to be seriously damaged by various random interference factors when operating in ordinary sea areas. Therefore, it becomes very necessary to apply the bioprotein material we designed to soft underwater robots. This self-healing material will repair the damaged parts during each collision, reduce the damage rate, and thereby extend the lifespan. Most importantly, this will cut the cost of salvage and maintenance.

In order to enhance the self-healing function of biomaterials, we designed and constructed a circular mRNA that can generate highly repetitive amino acid sequences rich in β-sheet. A single translation of the circular mRNA generates TRn5. In order to test the stability of proteins formed by different translation times of circular mRNA, we performed mathematical and biological simulations. In mathematical simulation, we described the stability of the protein through the hydrogen bond network formed between β-sheet.

First, we selected proteins translated from the circular mRNA 1, 2, 3, 4, 5, and 6 times respectively to study the stability of the intermolecular hydrogen bond networks. In the modeling process, we placed several points in a grid, with each point representing a β-sheet. In order to more realistically reflect the distribution of the β-sheet in the protein, we set the concentration of the points according to the relevant parameters of the grid. In order to enable different points to interact with each other, we set the connection rules according to the bond length of the hydrogen bond. If there is a connection, it means that the two points interact through hydrogen bonds. After the first hydrogen bond network was formed, we perturbed the system 20 times, mainly setting perturbations on the position of the points to simulate the reaction of the protein when the environment changes, and then regenerated a new hydrogen bond network. For each of the 20 hydrogen bond networks formed, we calculated the point-to-line ratio and used the average value as an evaluation index for the stability of the protein at each translation level.

The ratio of the number of points to the number of lines can reflect the stability of the network structure to a certain extent. When the ratio of points to lines is low, it means that the number of lines is relatively large to the number of points, and the network is more tightly connected. Networks with dense lines usually have higher stability and can more effectively resist external interference. In addition, such networks are not prone to deformation due to the dense distribution of their lines. When a line is broken, other surrounding lines can share its load and reduce the impact of the break on the overall structure. Therefore, in microscopic protein structures designed for self-healing functional materials, networks with dense lines can often recover faster and show stronger self-healing ability. The figure shows the changes and result analysis of the intermolecular hydrogen bond network formed by proteins under different translation times are shown in Fig. 15 and Fig. 16a.

Fig. 15 | Intermolecular hydrogen bond networks in different proteins

TRn with n=5, 10, 15, 20, 25 and 30 from top to bottom, left to right.

Fig. 16 | Stability analysis of different results

From Fig. 16a, we can see that as the number of TRn repetitions increases, the point-to-line ratio of the hydrogen bond network gradually decreases, indicating that the network transitions from a loose to a tight structure. Specifically, when the number of TRn repetitions increases from 5 to 10 times, the point-to-line ratio drops sharply from 1.2167 to 0.3921, which means that the density of the connections increase significantly, and the stability of the formed network structure converges rapidly. When the number of repetitions further increases to 15, 20, 25 and 30 times, the rate of decrease in the point-to-line ratio slows but continues to decline, falling to 0.3247, 0.3157, 0.2540 and 0.2109, respectively. This shows that a further increase in the number of repetitions will still improve the stability of the formed structure. In summary, proteins with a greater number of repetitions can form a more stable hydrogen bond network. This stable hydrogen bond network not only helps to maintain the integrity of the protein structure, but also improves the ability of the protein to bind to other molecules in biological processes, thereby promoting the realization of more complex biological functions. Therefore, the more times circular mRNA is translated, the higher the biological activity and functional stability of the protein.

Subsequently, we introduced several proteins with the same number of translations and placed them in a three-dimensional space to simulate the overall hydrogen bond network. We represented the three-dimensional space as a 'box' in the molecular dynamics simulation, abstracted the protein molecule into a sphere and set the size of the sphere according to the molecular weight, and abstracted the β-sheet as a point. Then we placed four spheres in the space and placed several evenly distributed points in the sphere according to the number of β-sheet. When placing the spheres, in order to improve the accuracy and rationality of the simulation, we limited the distance between different spheres to a certain range and constrained the distance from the sphere to the boundary of the space to prevent the protein from running out of the 'box' during the simulation. Relevant literature shows that the probability of hydrogen bonds between β-sheet in the same molecule is greater than the probability of bonds between different molecules. Therefore, we set different connection probabilities based on whether the two points belong to the same sphere and combined them with the length of the hydrogen bond to formulate the initial connection rules. When simulating the structure of the intermolecular hydrogen bond network, we did not consider the breaking of hydrogen bonds due to their high strength. However, in the simulation of the overall hydrogen bond network, we determined the probability of bond breaking based on the degree. The degree of a point is the number of lines that extend from the point. If a β-sheet forms more hydrogen bonds, the strength of the bond will be relatively weakened. After forming the innitial hydrogen bond network, we perturbed the system for 20 times. The results after each perturbation were based on the previous network. We determined bond formation and dissociation based on the degree of the two points, the spheres to which they belonged, and the distance between them. Finally, we calculated the point-to-line ratio for each hydrogen bond network formed, and took its average and variance as evaluation indicators for the stability of systems composed of certain proteins. The simulation display and result analysis are shown in Fig. 17 and Fig. 16b.

Fig. 17 | Simulation of hydrogen bond networks between TRn molecules

TRn with n=5, 10, 15, 20, 25 and 30 from top to bottom, left to right.

As can be seen from Fig. 16b, the point-to-line ratio of the hydrogen bond network between multiple protein molecules is generally lower than the ratio observed within a single protein, indicating that the hydrogen bond interactions between multiple proteins are more intensive and the connectivity of the formed network is stronger. Specifically, as the number of protein repetitions increases, the point-to-line ratio of the network continues to decrease, from 0.3461 to 0.0687, reflecting the increasing density of the connection line of the structure. When the circular mRNA is translated once to generate TRn5, the hydrogen bond network at this time has a higher point-to-line ratio, indicating relative sparsity and poor resistance to interference. As TRn repetitions increase, the point-to-line ratio of the hydrogen bond network gradually decreases. When the number of repetitions reaches 30, the point-to-line ratio is only 0.0687, indicating that the overall structure tends to be saturated. At this time, although there is still the possibility of new connection formation, most of the points in the structure interactions have maintained stable through sufficient connections, and the impact of further perturbations on the structure becomes smaller and smaller. In summary, the more TRn is repeated, the more stable the hydrogen bond network between protein molecules. This can make the interacting protein molecules more resistant to deformation when facing external perturbations, and have higher self-healing ability and adaptability.

We then used GROMACS to perform molecular dynamics simulations and calculate the energy in the system, hoping to verify the validity of our mathematical modeling through molecular dynamics methods. However, due to the large molecular weight of our system and limited computing resources, we only simulated TRn20 at most.

It should be noted that since GROMACS could not directly calculate the energy of hydrogen bonds in the system, we directly used the total energy of the system (including the energy contained in van der Waals forces and electrostatic forces) to reflect the stability of the network, because the total energy included the energy generated by hydrogen bonds, while the energy from sources other than hydrogen bonds was basically the same for each TRn. In addition, to eliminate the influence of protein molecule quantity on the results, we added the same qu of protein molecules to each system. These molecules were composed of TRn with different numbers of repetitions (for example, 5 TRn5 and 5 TRn10). When calculating the energy, we divided the total energy by the quantity of molecules and the number of repetitions of TRn, ultimately obtaining the energy corresponding to a unit TRn in the system. We used this energy as a measure of system stability, and the result is shown in Fig. 16c.

From the results, we can see that as the number of TRn repetitions increases, the energy of the system becomes lower and lower, so the protein becomes more stable, which is consistent with the results we obtained through mathematical modeling, and to a certain extent, it shows the accuracy of the mathematical modeling results. Due to the above conditions, we only use molecular dynamics as an auxiliary method to verify our model, because the results it obtains are correct in trend, but there may be errors in the numerical values. In the future, we will use other molecular dynamics methods and tools to separately calculate the hydrogen bond energy contained in the system to optimize it.

From the results, we can see that as the number of TRn repetitions increases, the energy of the system becomes lower and lower, so the protein becomes more stable, which is consistent with the results we obtained through mathematical modeling, and to a certain extent, it shows the accuracy of the mathematical modeling results.

Due to the above conditions, we only use molecular dynamics as an auxiliary method to verify our model, because the results it obtains are correct in trend, but there may be errors in the numerical values. In the future, we will use other molecular dynamics methods and tools to separately calculate the hydrogen bond energy contained in the system to optimize it.

L-DOPA is a crucial product of our project, derived from the conversion of tyrosine catalyzed by tyrosinase. The higher the efficiency of converting tyrosine to is L-DOPA, the better our protein performs. Nevertheless, L-DOPA is prone to further oxidation, which can lead to a decline in production. Therefore, we sought to identify an optimal reaction time that maximized both the yield of L-DOPA and the efficiency of the reaction. Through searching references, we have acquired the detailed reaction pathway for the formation of L-DOPA.

Fig. 18 | Pathway for the formation of dopa

Next, we conducted differential equation modeling for this reaction system. With the aid of the Michaelis-Menten equation, we can derive the following three equations, where [*] represents the concentration of a certain substance within the reaction system.

$$\tag{1} \frac{d\left[ Tyr \right]}{dt}=-\frac{v_{m_1}\left[ Tyr \right]} {k_{m_1}+\left[ Tyr \right]} \quad \quad $$ $$ \tag{2}\frac{d\left[ L-DOPA \right]}{dt}=\frac{v_{m_1}\left[ Tyr \right]}{k_{m_1}+\left[ Tyr \right]}-\frac{v_{m_2}\left[ L-DOPA \right]}{k_{m_2}+\left[ L-DOPA \right]}+\frac{v_{m_3}\left[ DQ \right]}{k_{m_3}+\left[ DQ \right]} $$ $$ \tag{3}\frac{d\left[ DQ \right]}{dt}=\frac{v_{m_2}\left[ L-DOPA \right]}{k_{m_2}+\left[ L-DOPA \right]}-\frac{v_{m_3}\left[ DQ \right]}{k_{m_3}+\left[ DQ \right]} $$

There are two types of parameters in these three equations, namely$K_{m_{\alpha}}$ and$V_{m_{\alpha}}$. Among them, $K_{m_{\alpha}}$ represents the reaction constant of the Michane of the response , and $ V_ {m _ {\ alpha}} $ represents the maximum reaction rate of the reaction, which is related to the enzyme concentration.

Based on the data from wet lab, we selected the optimal fitting parameter group by comparing several sets of parameters. We achieved this by calculating the sum of squares for the discrepancies between the numerical solutions obtained from the differential equations and the actual experimental measurements, the results of fit are shown in the following Fig. 19 (since the maximum absorption peaks of tyrosine and L-DOPA were very similar and insignificant, only the variation of L-DOPAquinone concentration was measured):

Fig. 19 | Fitting effects

In this case, the goodness of fit $R^2$ is 0.9962, indicating that it fits well.

Subsequently, we applied the fitting parameters into the model and set the initial Tyr concentration as 0.16 μmol/L. The concentration changes of each substance in the reaction system are shown in Fig. 20.

Fig. 20 | Changes in the concentration of substances in the reaction system

From the Fig. 20, we found that the production of L-DOPA in the reaction system showed a tendency of increasing and then decreasing, which reached the highest value at about 130 s, and then decreased gradually. Meanwhile, the production of L-DOPAquinone continued to increase, indicating that L-DOPA was gradually converted into L-DOPAquinone. Therefore, the optimum reaction time is approximately 130 seconds, when the L-DOPA production peaks.

To demonstrate that the reaction can perform stably in a general environment, we added a disturbance in each reaction channel. The disturbance represented the changes of temperature, humidity and other environmental factors in the interior. We took seven values of the deviation from standard environment as the disturbance frequency, such as -15%, -10%, -5%, 0, 5%, 10% and 15%, and then changed the initial tyrosine concentration. The variation of L-DOPA production under different conditions was shown in Fig. 21.

Fig. 21 | concentration of substances in the system after disturbing

As we can see in Fig. 21, under different disturbance, the trend of L-DOPA quantity is similar and the production fluctuates slightly, indicating that the reaction can be carried out stably in different environments. Therefore, it is concluded that our reaction system has greater capacity of environment adaptability and stability.

Table 4 | Symbol Description

When soft robots are working, even if they are applied with protective materials, regular maintenance is needed. And with the frequent increase in the number of uses and repairs, the equipment in the engineering practice will gradually show the phenomenon of "repair fatigue", that is, the marginal contribution of maintenance activities to improve equipment span and failure rate will become smaller, so regular replacement is also necessary. These replacements and maintenance will both incur certain costs.

The self-healing material we designed can save the subsequent maintenance cost of soft robots and their protective materials, which have high economic value. We established a model to quantify the economic benefits brought by our self-healing material, indicating that it has high practical value.

First, we calculated the cost of maintaining soft robots under normal circumstances.

Maintenance can also be divided into two main categories: corrective and preventive maintenance. Corrective maintenance (CM) is maintenance performed when the soft robot is injured. Preventive maintenance (PM) is maintenance performed while the system is running to keep the soft robot in a superior operating condition by inspecting, detecting, and preventing initial failures. Preventive maintenance is essential because some impairment is difficult to detect and repair during the course of a mission.

In general, robots may arrive at the task location at random times because in many cases it is difficult for the decision-makers to control or predict their arrival time, so we should take the time that robot travel to the task location into account. However, the duration of each task is under our control, so it would be interesting to estimate the possible task durations and provide optimal replacement strategies.

In our model, we used a random maintenance policy, which treats the time intervals between injuries to the soft robot and the time it takes for the robot to arrive at the workplace as random variable. We only considered the case where the soft robot suffers a non-fatal injury because we wanted to focus on the economic value of the self-healing ability of the self-healing material rather than its protective ability. We stipulated that at the end of the mission, the soft robot would be salvaged from the shore and underwent a corrective repair (CM) of the damage it sustained during the mission at a cost of $c_f$, while a preventive maintenance (PM) will be performed to replace the robot or its original parts at a cost of $c_p$.

We set that the time between each time the soft robot suffers injury that needs to be repaired as a random variable $X$ with density function $f(t \ )$ and distribution function $F (t\ )$. Therefore, the above process can be considered as a renewal process, and set its number of occurrences in (0,t] time be $N(t\ )$. We denote its update function as $M(t\ )$, representing the expected value of the number of occurrences in the time (0,t], i.e., $M(t\ )=E\left[ N(t\ ) \right] =\sum\nolimits_{n = 1}^{\infty}{F ^ {(n\ )}(t\ )}$, and let $m(t\ )=\frac{dM (t\ )}{dt}$ be the renewal density

In addition, we assumed that the time required for the soft robot to arrive at the workplace each time is $T$, which is also a random variable with a distribution function $G(t\ )$, a density function$ g(t\ )$. And the duration of a task is $t_x$.

According to the renewal reward theorem, after a long time, the expected cost rate is:: $$ \tag{1} C (t_x)=\frac{c_f\int_0^{\infty}{M (t + t_x)dG (t\ )+c_p}}{\int_0^{\infty}{ (t + t_x)dG (t\ )}} $$

We aim to find the optimal working time $t_x$ to minimize the cost rate function $C(t\ )$.

By Differentiating $C(t\ )$ with respect to $t_x$ and setting it equal to zero, we obtain the the following equation: $$ \tag{2} \int_0^{\infty}{m (t + t_x\ )dG (t\ )\cdot \int_{0}^{\infty}{ (t + t_x\ )dG(t\ )-\int_0^{\infty}{M (t + t_x\ )dG(t\ )}}}=\frac{c_p}{c_f} $$

This is a necessary but not sufficient condition for $t_x$ to be able to minimise the cost, so we had to further analyse the values of $t_x $. We denoted the left-hand side of the above equation as $Z(t_x)$, i.e.

$$ \tag{3} Z\left( t_x \right) =\int_{0}^{\infty}{m(t + t_x)dG(t\ )\cdot \int_0^{\infty}{(t + t_x)dG(t\ )-\int_0^{\infty}{M(t + t_x)dG(t\ )}}} $$

It is easy to know that when $m(t\ )$ is nondecreasing function, $Z (t_x)$ is also a nondecreasing function, i.e., $Z (t_x)$ increases from $\int_0^{\infty}{m(t\ )\left[ \theta g(t\ )-\overline{G}(t\ ) \right] dt}$ to $Z\left( \infty \right) $

Therefore, we had the following optimum policies:

1.If $\int_0^{\infty}{m(t\ )\left[ \theta g(t\ )-\overline{G}(t\ ) \right] dt}\ge \frac{c_p}{c_f}$, since $Z(t_x)$ is a increasing function, $Z(t_x)$ is always greater than $\frac{c_p}{c_f}$, and $C(t_x)$ is minimized when $t_x$ is equal to 0. This implies that the cost of the soft robot increases with the length of each working time, and if we want to save the cost, we need to make the working span of the soft robot as short as possible.

2.$\int_0^{\infty}{m(t\ )\left[ \theta g(t\ )-\overline{G}(t\ ) \right] dt} \lt \frac{c_p}{c_f} \lt z \left( \infty \right)$, solving Equation (2), we can get $t_x$ to minimize cost function $C(t\ )$ and the minimum cost is

$$ \tag{4} C(t_{x}^{\ *})=c_f\int_0^{\infty}{m(t + t_{ x } ^ {\ *})}dG(t\ ) $$

3.If $\frac{c_p}{c_f}\geqslant z\left( \infty \right)$, then $Z(t_x)$ is always less than $\frac{c_p}{c_f}$, that is to say, the longer the work span is, the less costly it will be, but in practice, we also need to consider that long hours of work lead to an increased probability of fatal injury.

Poisson process is a common distribution in nature, especially with regard to counting process, usually the counting process that satisfies the following conditions may be a Poisson process:

1. A single event is independent of whether it occurs or not, and of the probability of its occurrence;

2. It is known that the average number of times an event occurs (incidence) occurs in a given interval (time/space) is finite.

We can assume that the number of damage suffered by the submarine robot while working meets the above conditions, so we can consider it as a Poisson process. From the knowledge of probability theory, $f(t)$ is the density function of the exponential distribution, i.e., $f( x \ ) =\lambda e^{-\lambda t}\,\, \left( t\geqslant 0 \right) $, and also $M( t \ ) =\lambda t$, $m( t \ ) =\lambda$.

We substitute the above results into Equation (3) to get

$$ \tag{5} Z(0)=\int_0^{\infty}{\lambda \left[ \theta g(t\ )-\overline{G}(t\ ) \right]}dt $$

When $G(t\ )={-\sigma t}\,\, (0 \lt \sigma \lt \infty )$, we can get $Z(t_x)=\lambda t_x$, which satisfies $Z(0)\lt \frac{c_p}{c_f} \lt Z\left( \infty \right)$. Thus we can solve Equation (2) as:

$$ \tag{6} t_x=\frac{c_p}{\lambda c_f} $$

The minimized cost rate is:

$$ \tag{7} \lambda c\,\,_f+\frac{\sigma \lambda c_pc_f}{\lambda c_f+\sigma c_p} $$

If we assumed that in an average of two hours at the bottom of the ocean, the soft robot would be injured and need to be repaired, then $\lambda =0.5$. Assuming that each corrective maintenance will cost 400 dollars and preventive maintenance will cost 600 dollars , $\sigma =1$. Based on the above derivation, the optimal strategy would be to have a duration of 3 hours per mission, with an average minimum cost of 350 dollars per hour and the minimum cost rate of 1050 dollars per mission.

self-healing materials do not need to be repaired after being damaged, which means that under the conditions we set, self-healing materials can save 450 dollars in costs during each mission. Although self-healing materials may be more expensive than ordinary materials, after long-term use, the cost savings of self-healing materials will make up for this gap. And the biggest advantage of self-healing materials is that we do not need to salvage the robot for maintenance during the mission, which can extend the duration of each mission, allowing the robot to perform longer missions and bring more benefits.

Of course, due to various restrictions, we assumed some data. In actual production, using collected real data can make the model closer to reality.